Altering social network position prevents social isolation

Does allocating cooperative individuals to high-degree nodes prevent isolation?

round = NULL
gameID = NULL 
nodes = NULL
meanDegree = NULL
meanDegreeCoop = NULL #mean degree among those who cooperated
meanDegreeDefect = NULL #mean degree among those who defected
meanDist = NULL #mean shortest path length divided by the maximum possible path length (Nvertices - 1)
maxDist = NULL #max shortest path length divided by the maximum possible path length (Nvertices - 1)
coefGlobal = NULL #global clustering coefficient
percentCoop = NULL #proportion of those who cooperated
percentCoopPopular = NULL #proportion of those who cooperated, among those with degree >=10
coopEnv = NULL 
assort = NULL

for(i in 1:length(g.isolation)){
    round[i] =  mean(V(g.isolation[[i]])$round)
    gameID[i] = V(g.isolation[[i]])$gameID[1]
    nodes[i] = vcount(g.isolation[[i]])
    meanDegree[i] = mean(igraph::degree(g.isolation[[i]]),na.rm=TRUE)
    meanDegreeCoop[i] = mean(igraph::degree(g.isolation[[i]])[V(g.isolation[[i]])$behavior=="C"],na.rm=TRUE)
    meanDegreeDefect[i] = mean(igraph::degree(g.isolation[[i]])[V(g.isolation[[i]])$behavior=="D"],na.rm=TRUE)
    meanDist[i] = mean(distances(g.isolation[[i]])[distances(g.isolation[[i]])!=0 & distances(g.isolation[[i]])!=Inf],na.rm=TRUE) / length(V(g.isolation[[i]])-1)
    maxDist[i] = max(distances(g.isolation[[i]])[distances(g.isolation[[i]])!=0 & distances(g.isolation[[i]])!=Inf],na.rm=TRUE) / length(V(g.isolation[[i]])-1)
    coefGlobal[i] = transitivity(g.isolation[[i]], type="global")
    percentCoop[i] = prop.table(table(V(g.isolation[[i]])$behavior))["C"]
    percentCoopPopular[i] = prop.table(table(V(g.isolation[[i]])$behavior[igraph::degree(g.isolation[[i]])>=6]))["C"]
    coopEnv[i] = sum(ifelse(V(g.isolation[[i]])$behavior=="C",1,0)*igraph::degree(g.isolation[[i]])/sum(igraph::degree(g.isolation[[i]])),na.rm=TRUE)
    assort[i] = assortativity(g.isolation[[i]], V(g.isolation[[i]])$behavior == "C")
}

smallWorld = data.frame(
  round = round,
  gameID = gameID,
  nodes = nodes,
  meanDegree = meanDegree,
  meanDegreeCoop = meanDegreeCoop,
  meanDegreeDefect = meanDegreeDefect,
  diffDegreeCD = meanDegreeCoop - meanDegreeDefect,
  meanDist = meanDist,
  maxDist = maxDist,
  coefGlobal = coefGlobal,
  percentCoop = percentCoop,
  percentCoopPopular = percentCoopPopular,
  coopEnv = coopEnv,
  assort = assort
)


round = NULL
gameID = NULL
visibleWealth = NULL
n=1
for(i in unique(cdata$round)){
  for(m in unique(cdata$gameID)){
    round[n] = as.numeric(i)-1
    gameID[n] = m
    visibleWealth[n] = ifelse(cdata[cdata$round==i & cdata$gameID==m,]$showScore[1]=="true",1,0)
    n=n+1
  }
}

smallWorld.2 = data.frame(
  round = round,
  gameID = gameID,
  visibleWealth = visibleWealth
)

smallWorld = merge(smallWorld, smallWorld.2, by=c("gameID","round"), all.x=TRUE)

smallWorld.initial = subset(smallWorld, smallWorld$round==0)

#percent of people that end up being isolated for each gameID
percentIsolated = NULL
percentIsolated = aggregate(isolated ~ gameID, data=sample.wide, FUN=max, na.rm=TRUE)
smallWorld.initial = merge(smallWorld.initial[c("gameID","meanDegree","meanDegreeCoop","meanDegreeDefect","diffDegreeCD","meanDist","coefGlobal","visibleWealth","percentCoop","percentCoopPopular","coopEnv","assort")],percentIsolated, by="gameID", all.x=TRUE)

percentIsolated = NULL
percentIsolated = aggregate(isolated ~ gameID, data=sample.wide, FUN=mean, na.rm=TRUE)
percentIsolated = percentIsolated %>% rename(
  percentIsolated = isolated
)
smallWorld.initial = merge(smallWorld.initial,percentIsolated, by="gameID", all.x=TRUE)




reg = glm(isolated ~ coopEnv, 
      data = smallWorld.initial,
      family = binomial(link = "logit"))
summary(reg)

## 
## Call:
## glm(formula = isolated ~ coopEnv, family = binomial(link = "logit"), 
##     data = smallWorld.initial)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.6489  -0.4914  -0.4395  -0.3854   2.3384  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -0.8002     1.6395  -0.488    0.626
## coopEnv      -2.1872     2.5728  -0.850    0.395
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 52.013  on 79  degrees of freedom
## Residual deviance: 51.282  on 78  degrees of freedom
## AIC: 55.282
## 
## Number of Fisher Scoring iterations: 5

reg = glm(isolated ~ coopEnv*percentCoop, 
      data = smallWorld.initial,
      family = binomial(link = "logit"))
summary(reg)

## 
## Call:
## glm(formula = isolated ~ coopEnv * percentCoop, family = binomial(link = "logit"), 
##     data = smallWorld.initial)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.7581  -0.5129  -0.3847  -0.2979   2.5705  
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)
## (Intercept)           -6.739      8.192  -0.823    0.411
## coopEnv               -8.483     16.193  -0.524    0.600
## percentCoop           20.749     16.034   1.294    0.196
## coopEnv:percentCoop   -9.081     19.775  -0.459    0.646
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 52.013  on 79  degrees of freedom
## Residual deviance: 48.836  on 76  degrees of freedom
## AIC: 56.836
## 
## Number of Fisher Scoring iterations: 5

reg = glm(isolated ~ assort, 
      data = smallWorld.initial,
      family = binomial(link = "logit"))
summary(reg)

## 
## Call:
## glm(formula = isolated ~ assort, family = binomial(link = "logit"), 
##     data = smallWorld.initial)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.8087  -0.4639  -0.3327  -0.2150   2.6431  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -3.3112     0.7398  -4.476 7.62e-06 ***
## assort      -10.6817     4.6011  -2.322   0.0203 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 40.020  on 64  degrees of freedom
## Residual deviance: 33.269  on 63  degrees of freedom
##   (15 observations deleted due to missingness)
## AIC: 37.269
## 
## Number of Fisher Scoring iterations: 6

reg = glm(isolated ~ assort*percentCoop, 
      data = smallWorld.initial,
      family = binomial(link = "logit"))
summary(reg)

## 
## Call:
## glm(formula = isolated ~ assort * percentCoop, family = binomial(link = "logit"), 
##     data = smallWorld.initial)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.7953  -0.4511  -0.3392  -0.2241   2.6261  
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)
## (Intercept)         -3.4008     3.7606  -0.904    0.366
## assort             -14.8971    21.6781  -0.687    0.492
## percentCoop          0.2359     5.5394   0.043    0.966
## assort:percentCoop   7.1996    34.4835   0.209    0.835
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 40.020  on 64  degrees of freedom
## Residual deviance: 33.194  on 61  degrees of freedom
##   (15 observations deleted due to missingness)
## AIC: 41.194
## 
## Number of Fisher Scoring iterations: 6

reg = glm(percentIsolated ~ coopEnv, 
      data = smallWorld.initial, 
      family = gaussian(link = "identity"))
summary(reg)

## 
## Call:
## glm(formula = percentIsolated ~ coopEnv, family = gaussian(link = "identity"), 
##     data = smallWorld.initial)
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.013339  -0.008598  -0.006745  -0.004587   0.117036  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.01846    0.01166   1.583    0.118
## coopEnv     -0.01722    0.01732  -0.994    0.323
## 
## (Dispersion parameter for gaussian family taken to be 0.0005112187)
## 
##     Null deviance: 0.040380  on 79  degrees of freedom
## Residual deviance: 0.039875  on 78  degrees of freedom
## AIC: -375.29
## 
## Number of Fisher Scoring iterations: 2

reg = glm(percentIsolated ~ coopEnv*percentCoop, 
      data = smallWorld.initial, 
      family = gaussian(link = "identity"))
summary(reg)

## 
## Call:
## glm(formula = percentIsolated ~ coopEnv * percentCoop, family = gaussian(link = "identity"), 
##     data = smallWorld.initial)
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.014683  -0.009669  -0.006671  -0.003141   0.115792  
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)
## (Intercept)         -0.01736    0.05317  -0.326    0.745
## coopEnv             -0.03845    0.10335  -0.372    0.711
## percentCoop          0.11389    0.10197   1.117    0.268
## coopEnv:percentCoop -0.06003    0.12200  -0.492    0.624
## 
## (Dispersion parameter for gaussian family taken to be 0.0005144455)
## 
##     Null deviance: 0.040380  on 79  degrees of freedom
## Residual deviance: 0.039098  on 76  degrees of freedom
## AIC: -372.87
## 
## Number of Fisher Scoring iterations: 2

reg = glm(percentIsolated ~ assort, 
      data = smallWorld.initial, 
      family = gaussian(link = "identity"))
summary(reg)

## 
## Call:
## glm(formula = percentIsolated ~ assort, family = gaussian(link = "identity"), 
##     data = smallWorld.initial)
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.017732  -0.010336  -0.006163  -0.000788   0.119912  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.003475   0.003089   1.125   0.2649  
## assort      -0.064495   0.026242  -2.458   0.0167 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.0004921539)
## 
##     Null deviance: 0.033978  on 64  degrees of freedom
## Residual deviance: 0.031006  on 63  degrees of freedom
##   (15 observations deleted due to missingness)
## AIC: -306.66
## 
## Number of Fisher Scoring iterations: 2

reg = glm(percentIsolated ~ assort*percentCoop, 
      data = smallWorld.initial, 
      family = gaussian(link = "identity"))
summary(reg)

## 
## Call:
## glm(formula = percentIsolated ~ assort * percentCoop, family = gaussian(link = "identity"), 
##     data = smallWorld.initial)
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.017938  -0.009107  -0.005072  -0.002206   0.119805  
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)  
## (Intercept)         0.001917   0.016882   0.114   0.9100  
## assort             -0.232569   0.120569  -1.929   0.0584 .
## percentCoop         0.002524   0.024046   0.105   0.9167  
## assort:percentCoop  0.257715   0.178889   1.441   0.1548  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.0004884538)
## 
##     Null deviance: 0.033978  on 64  degrees of freedom
## Residual deviance: 0.029796  on 61  degrees of freedom
##   (15 observations deleted due to missingness)
## AIC: -305.24
## 
## Number of Fisher Scoring iterations: 2

# #MSM
# # estimation of denominator of ip weights
# #cooperatorの数以外は全てRandom（Erdos-Renyi random graphにより決まるため）なので、coopEnvはpercentCoopで調整すれば十分か
# den.fit.obj <- lm(coopEnv ~ percentCoop + I(percentCoop^2), data = smallWorld.initial)
# p.den <- predict(den.fit.obj, type = "response")
# dens.den <- dnorm(smallWorld.initial$coopEnv, p.den, summary(den.fit.obj)$sigma)
# 
# # estimation of numerator of ip weights
# num.fit.obj <- lm(coopEnv ~ 1, data = smallWorld.initial)
# p.num <- predict(num.fit.obj, type = "response")
# dens.num <- dnorm(smallWorld.initial$coopEnv, p.num, summary(num.fit.obj)$sigma)
# 
# smallWorld.initial$sw.a = dens.num/dens.den
# summary(smallWorld.initial$sw.a)
# 
# msm.sw.cont <-geepack::geeglm(isolated ~ coopEnv, 
#                  data=smallWorld.initial, weights=sw.a, family = binomial(link="logit"), id=gameID, corstr="independence")
# summary(msm.sw.cont)
# 
# msm.sw.cont <-geepack::geeglm(percentIsolated ~ coopEnv, 
#                  data=smallWorld.initial, weights=sw.a, family = gaussian(link="identity"), id=gameID, corstr="independence")
# summary(msm.sw.cont)


ggplot(data = smallWorld.initial, aes(x = coopEnv, y = percentIsolated)) + 
  geom_point() +
  ggtitle("Percent isolated") +
  scale_x_continuous("Cooporative environment") +  
  scale_y_continuous("Percent of individuals isolated") + 
  scale_fill_continuous(type = "viridis")

ggplot(data = smallWorld.initial, aes(x = assort, y = percentIsolated)) + 
  geom_point() +
  ggtitle("Percent isolated") +
  scale_x_continuous("Assortativity coefficient of initial C/D") +  
  scale_y_continuous("Percent of individuals isolated") + 
  scale_fill_continuous(type = "viridis")

## Warning: Removed 15 rows containing missing values (`geom_point()`).

Simulation 1

All graphs start as Erdos-Renyi random newtorks

Cooperators have low degrees

#Define degrees of isolation
isolationDegree = 0

#number of iterations per arm
iterations = 300

modelForPrediction = "random forest" #"linear" or "random forest"

# List of manipulating parameters of experiments
#L : number of rounds
#V : Visible or not
#A : Income of a rich-group subject
#B : Income of a poor-group subject
#R : Probability to be assigned to a rich group
#I : Number of the same-parameter trial

R = 0.5
I = 0
L = 10


trends.df = data.frame()
  

  
  
for(A in c(1150,700,500)){
  for(V in c(0,1)){
      
      V = V
      A = A
      if(A==1150){B = 200} #high inequality
      if(A==700){B = 300} #low inequality
      if(A==500){B = 500} #no inequality
  
      if(modelForPrediction=="random forest"){
        source(paste(rootdir,"R/models.R",sep="/"))
        if(V==0){
          model1<-model1.invisible(redo=FALSE)
          model2<-model2.invisible(redo=FALSE)
          model3<-model3(redo=FALSE)
          }
        if(V==1){
          model1<-model1.visible(redo=FALSE)
          model2<-model2.visible(redo=FALSE)
          model3<-model3(redo=FALSE)
          }
      }
      
      df.netIntLowDegree = data.frame(
        coopFrac = NULL,
        avgCoop = NULL, 
        avgCoopFinal = NULL, 
        percentIsolation = NULL,
        isolation = NULL,
        percentIsolationC = NULL,
        percentIsolationD = NULL,
        nCommunities = NULL,
        communitySize = NULL,
        assortativityInitial = NULL,
        assortativityFinal = NULL,
        conversionRate = NULL,
        conversionToD = NULL, 
        conversionToC = NULL, 
        transitivity = NULL, 
        degree = NULL, 
        degreeC = NULL, 
        degreeD = NULL,
        meanConversionToD = NULL, 
        meanConversionToC = NULL, 
        degreeLost = NULL,
        degreeLostC = NULL,
        degreeLostD = NULL
      )
  
      #Here, factionCoop=0 will be the control: no rearranging of nodes will take place
      for(frac in c(0,0.25,0.5,0.75,1)){
        #nodes in the top fractionCoop degrees will automatically be a cooperator
        fractionCoop = frac
        
          coopFrac = NULL
          avgCoop = NULL
          homophilyC = NULL
          homophilyD = NULL
          heterophily = NULL
          avgCoopFinal = NULL
          percentIsolation = NULL
          isolationPersonRounds = NULL
          isolation = NULL
          percentIsolationC = NULL
          percentIsolationD = NULL
          nCommunities = NULL
          communitySize = NULL
          assortativityInitial = NULL
          assortativityFinal = NULL
          conversionRate = NULL
          conversionToD = NULL 
          conversionToC = NULL
          transitivity = NULL
          degree = NULL
          degreeC = NULL
          degreeD = NULL
          meanConversionToD = NULL 
          meanConversionToC = NULL
          degreeLost = NULL
          degreeLostC = NULL
          degreeLostD = NULL
          avg_wealth = NULL
          gini = NULL
          for(m in c(1:iterations)){
            # Section 1. NOTES, packages, and Parameters
            #Importing library
            library(igraph) # for network graphing
            library(reldist) # for gini calculatio
            library(boot) # for inv.logit calculation
            #Two prefixed functions
            #rank
            rank1 = function(x) {rank(x,na.last=NA,ties.method="average")[1]} #a smaller value has a smaller rank.
            #gini mean difference (a.k.a. mean difference: please refer to https://stat.ethz.ch/pipermail/r-help/2003-April/032782.html)
            gmd = function(x) {
             x1 = na.omit(x)
             n = length(x1)
            tmp = 0
             for (i in 1:n) {
             for (j in 1:n) {
             tmp <- tmp + abs(x1[i]-x1[j])
             }
             }
            answer = tmp/(n*n)
            if(length(x)!=0){return(answer)}
            if(length(x)==0){return(NA)}
            }
        
            
            
            
            # List of fixed parameters of experiments (assumptions)
            #Rewiring rate = 0.3
            #GINI coefficient (can be known by A or B)
            GINI = 0*as.numeric(A==500) + 0.2*as.numeric(A %in% c(700,850)) + 0.4*as.numeric(A ==1150)
            #Collecting data frame (final output data frame)
            result = data.frame(round=0:L,n_par=NA,n_A=NA,avg_coop=NA,avg_degree=NA,avg_wealth=NA,gini=NA,gmd=NA,avg_coop_A=NA,avg_degree_A=NA,avg_wealth_A=NA,gini_A=NA,gmd_A=NA,avg_coop_B=NA,avg_degree_B=NA,avg_wealth_B=NA,gini_B=NA,gmd_B=NA,isolation=NA,percentIsolation=NA,meanConversionToD=NA,meanConversionToC=NA,degreeLost=NA,degreeLostC=NA,degreeLostD=NA)
            #_A is for a richer group and _B is for a poorer group
            
            
            #####################################################
            # Section 1.5: Practice rounds 1 to 2, to determine C/D in round 1
            
            N = 8 # median of the number of participants over rounds.
            node_rp0 = data.frame(ego_id=1:N, round=0)
            node_import = node_rp0
            
            for (k in 1:2){
              node_rX = node_import #Importing data
              node_rX$round = node_rX$round + 1
              node_rX[is.na(node_rX$prev_degree)==1,"prev_degree"] = 0
              node_rX[is.na(node_rX$prev_local_rate_coop)==1,"prev_local_rate_coop"] = 0
              #Only this calculation needs to change from Round 1
              if (k==1) {
                node_rX$prob_coop = inv.logit(1.099471) 
              } else {
                node_rX$prob_coop = inv.logit((-0.02339288) + (1.46068980)*as.numeric(node_rX$prev_coop==1))
              } 
              
              node_rX$coop = apply(data.frame(node_rX$prob_coop),1,function(x) {sample(1:0,1,prob=c(x,(1-x)))})
              
              node_rX$prev_coop = node_rX$coop
                
              assign(paste("coop_rp",k, sep=""),node_rX$coop)
              
              #For the loop
              node_import = node_rX
            }
            
            #cooperation rate in the practice rounds
            coop_rp1 <- rbinom(n = N, size = 1, prob = 0.6)
            coop_rp2 <- abs(coop_rp1 - rbinom(n = N, size = 1, prob = 0.1))
            coop_rp = apply(cbind(coop_rp1,coop_rp2),1,mean)
      
            #####################################################
            # Section 2: Round 0 (Agents and environments)
            #Node data generation
            N = 8 # median of the number of participants over rounds.
            node_r0 = data.frame(ego_id=1:N, round=0)
            node_r0$coop_rp = ifelse(coop_rp==1,"C","D")
            node_r0$group = sample(c("rich","poor"),N,replace=TRUE,prob=c(R,1-R)) #R is defined as the probability to be assigned to the rich group
            node_r0$initial_wealth = ifelse(node_r0$group=="rich",A,B)
            #Link data generation
            ego_list = NULL
            for (i in 1:N) { ego_list = c(ego_list,rep(i,N)) }
            link_r0 = data.frame(ego_id=ego_list,alt_id=rep(1:N,N))
            link_r0 = link_r0[(link_r0$ego_id < link_r0$alt_id),] #The link was bidirectional, and thus the half and self are omitted.
            link_r0$connected = sample(0:1,dim(link_r0)[1],replace=TRUE,prob=c(0.7,0.3)) #Initial rewiring rate is fixed, 0.3
            link_r0c_ego = link_r0[link_r0$connected==1,]
            link_r0c_alt = link_r0[link_r0$connected==1,]
            colnames(link_r0c_alt) = c("alt_id","ego_id","connected")
            link_r0c = rbind(link_r0c_ego,link_r0c_alt) #this is bidirectional (double counted) for connected ties.
            link_r0c = link_r0c[order(link_r0c$ego_id),]
            link_r0c$alternumber = NA #putting the number for each alter in the same ego
            link_r0c[1,]$alternumber = 1
            for (i in 1:(dim(link_r0c)[1]-1))
              {if (link_r0c[i,]$ego_id == link_r0c[i+1,]$ego_id)
                {link_r0c[i+1,]$alternumber = link_r0c[i,]$alternumber + 1}
              else
                {link_r0c[i+1,]$alternumber = 1}
              #print(i)
              }
            link_r0c2 = reshape(link_r0c, direction = "wide", idvar=c("ego_id","connected"), timevar="alternumber")
            if(class(link_r0c2[,colnames(link_r0c2)[substr(colnames(link_r0c2),1,6) == "alt_id"]]) == "data.frame"){
              link_r0c2$initial_degree = apply(link_r0c2[,colnames(link_r0c2)[substr(colnames(link_r0c2),1,6) == "alt_id"]],1,function(x){length(na.omit(x))})
              } #Degree of each ego
            if(class(link_r0c2[,colnames(link_r0c2)[substr(colnames(link_r0c2),1,6) == "alt_id"]]) == "integer"){
              link_r0c2$initial_degree = sapply(link_r0c2[,colnames(link_r0c2)[substr(colnames(link_r0c2),1,6) == "alt_id"]],function(x){length(na.omit(x))})
              } #Degree of each ego
            #Reflect the degree and initial local gini coefficient into the node data
            node_r0 = merge(x=node_r0,y=link_r0c2,all.x=TRUE,all.y=FALSE,by="ego_id")
            node_r0[is.na(node_r0$initial_degree)==1,"initial_degree"] = 0
            node_r0$initial_avg_env_wealth = NA
            node_r0$initial_local_gini = NA #local gini coefficient of the ego and connecting alters
            node_r0$initial_rel_rank = NA #local rank of ego among the ego and connecting alters (divided by the number of the go and connecting alters)
            for (i in 1:(dim(node_r0)[1])){
              node_r0[i,]$initial_avg_env_wealth = mean(na.omit(node_r0[node_r0$ego_id %in%
              node_r0[i,colnames(node_r0)[substr(colnames(node_r0),1,6) %in% c("ego_id","alt_id")]],"initial_wealth"]))
              node_r0[i,]$initial_local_gini = gini(na.omit(node_r0[node_r0$ego_id %in% node_r0[i,colnames(node_r0)[substr(colnames(node_r0),1,6)
              %in% c("ego_id","alt_id")]],"initial_wealth"]))
              node_r0[i,]$initial_rel_rank = rank1(na.omit(node_r0[node_r0$ego_id %in% node_r0[i,colnames(node_r0)[substr(colnames(node_r0),1,6)
              %in% c("ego_id","alt_id")]],"initial_wealth"]))/length(na.omit(node_r0[node_r0$ego_id %in%
              node_r0[i,colnames(node_r0)[substr(colnames(node_r0),1,6) %in% c("ego_id","alt_id")]],"initial_wealth"]))
              }
            #Finalization of round 0 and Visualization
            #plot(graph.data.frame(link_r0[link_r0$connected==1,],directed=F)) #plot.igraph
            node_r0$everIsolated = 0
            node_r0$maxDegreeLost = NA
            result[result$round==0,2:25] = c(length(node_r0$ego_id),length(node_r0[node_r0$group=="rich",]$ego_id),NA,mean(node_r0$initial_degree),mean(node_r0$initial_wealth),gini(node_r0$initial_wealth),gmd(node_r0$initial_wealth),NA,mean(node_r0[node_r0$group=="rich",]$initial_degree),mean(node_r0[node_r0$group=="rich",]$initial_wealth),gini(node_r0[node_r0$group=="rich",]$initial_wealth),gmd(node_r0[node_r0$group=="rich",]$initial_wealth),NA,mean(node_r0[node_r0$group=="poor",]$initial_degree),mean(node_r0[node_r0$group=="poor",]$initial_wealth),gini(node_r0[node_r0$group=="poor",]$initial_wealth),gmd(node_r0[node_r0$group=="poor",]$initial_wealth),
                                             as.numeric(ifelse(is.na(table(node_r0$initial_degree<=isolationDegree)["TRUE"]),0,1)),
                                             as.numeric(sum(node_r0$everIsolated)/length(node_r0$ego_id)),
                                             NA,
                                             NA,
                                             NA,NA,NA
                                             )
            
            #For the loop at the next round (for round 1, the initial one is the same as the previous [1 prior] one)
            node_import = node_r0
            node_import$initial_coop = NA
            node_import$prev_coop = NA
            node_import$prev_wealth = node_import$initial_wealth
            node_import$prev_degree = node_import$initial_degree
            node_import$prev_avg_env_wealth = node_import$initial_avg_env_wealth
            node_import$prev_local_gini = node_import$initial_local_gini
            node_import$prev_rel_rank = node_import$initial_rel_rank
            node_import$prev_local_rate_coop = NA
            link_import = link_r0
            
            
            #####################################################
            # Section 3: Rounds 1 to 10 or more (behaviors in simulation: the equation of cooperation is different at round 1 because of no history)
            #3-1: Cooperation phase
            for (k in 1:L)
            {
              node_rX = node_import #Importing data
              node_rX$round = node_rX$round + 1
              node_rX[is.na(node_rX$prev_degree)==1,"prev_degree"] = 0
              node_rX[is.na(node_rX$prev_local_rate_coop)==1,"prev_local_rate_coop"] = 0
              #Only this calculation needs to change from Round 1
              if(modelForPrediction=="linear"){
                if (k==1) {
                  node_rX$prob_coop = as.numeric(V==0)*inv.logit((-1.816665) + (2.086067)*coop_rp1 + (1.800153)*coop_rp2) + as.numeric(V==1)*inv.logit((-2.031577) + (2.427157)*coop_rp1 + (1.684193)*coop_rp2 + (-1.528851)*GINI)
                  } else {
                  node_rX$prob_coop = as.numeric(V==0 & node_rX$prev_coop==0)*inv.logit(-1.039916) + as.numeric(V==0 & node_rX$prev_coop==1)*inv.logit(2.062023) + as.numeric(V==1 & node_rX$prev_coop==0)*inv.logit((-0.2574838)*as.numeric(node_rX$prev_avg_env_wealth - node_rX$prev_wealth > 0) + (-1.214198)*GINI + (2.508148)*GINI*as.numeric(node_rX$prev_avg_env_wealth - node_rX$prev_wealth > 0) + (-0.9749075)) + as.numeric(V==1 & node_rX$prev_coop==1)*inv.logit((- 0.6197254)*as.numeric(node_rX$prev_avg_env_wealth - node_rX$prev_wealth > 0) + (-0.7480261)*GINI + (1.169674)*GINI*as.numeric(node_rX$prev_avg_env_wealth - node_rX$prev_wealth > 0) + (1.356784))
                  } 
              }
              if(modelForPrediction=="random forest"){
                if (k==1) {
                    if(V==1){node_rX$prob_coop = predict(model1,
                                                         newdata=
                                                           data.frame(
                                                             behavior.p1 = coop_rp1,
                                                             behavior.p2 = coop_rp2,
                                                             gini = GINI
                                                           ),
                                                         type = "prob"
                                                         )[[1]]$C}
                    else if(V==0){node_rX$prob_coop = predict(model1,
                                                         newdata=
                                                           data.frame(
                                                             behavior.p1 = coop_rp1,
                                                             behavior.p2 = coop_rp2
                                                           ),
                                                         type = "prob"
                                                         )[[1]]$C}
                  } else {
                    if(V==1){node_rX$prob_coop = predict(model2,
                                                         newdata=
                                                           data.frame(
                                                             prevCoop = node_rX$prev_coop,
                                                             gini = GINI,
                                                             alterPrevWealth = node_rX$prev_avg_env_wealth,
                                                             egoPrevWealth = node_rX$prev_wealth
                                                           ),
                                                         type = "prob"
                                                         )[[1]]$C}
                    else if(V==0){node_rX$prob_coop = predict(model2,
                                                         newdata=
                                                           data.frame(
                                                             prevCoop = node_rX$prev_coop,
                                                             alterPrevWealth = node_rX$prev_avg_env_wealth,
                                                             egoPrevWealth = node_rX$prev_wealth
                                                           ),
                                                         type = "prob"
                                                         )[[1]]$C}
                  }
              }
              #####rearrange node degrees before round 1 depending on cooperation in practice rounds!
              if(k==1){
                if(fractionCoop==0){
                  node_rX$prob_coop
                  node_rX$coop = apply(data.frame(node_rX$prob_coop),1,function(x) {sample(1:0,1,prob=c(x,(1-x)))})
                  coop_rp_init = coop_rp
                  }
                if(fractionCoop>0){
                  prob_coop_df = NULL
                  nodesCoop = NULL
                  #nodesCoop = node_rX$prev_degree<=quantile(node_rX$prev_degree,fractionCoop) #assign low-degree nodes to cooperators
                  #assign defectors to designated nodes
                  nodesCoop <- rep(FALSE, N)   
                  nodesCoop[order(node_rX$prev_degree)[(floor(N*(fractionCoop))-floor(N*0.25)+1):(floor(N*(fractionCoop)))]] <- TRUE #set true the nodes of assignment; we will select floor(N*0.25) nodes
                  prob_coop_df = 
                    data.frame(
                      prob_coop = rev(node_rX$prob_coop[order(coop_rp)]),
                      node_number = c(which(!nodesCoop),which(nodesCoop))
                    )
                  node_rX$prob_coop = prob_coop_df[order(prob_coop_df$node_number),]$prob_coop
                  #coop_rp of the rearranged nodes
                  coop_rp_init = rev(coop_rp[order(coop_rp)])[order(prob_coop_df$node_number)]
                  node_rX$coop = apply(data.frame(node_rX$prob_coop),1,function(x) {sample(1:0,1,prob=c(x,(1-x)))})
                }
              } else {
                node_rX$coop = apply(data.frame(node_rX$prob_coop),1,function(x) {sample(1:0,1,prob=c(x,(1-x)))})
              }
              if (k==1) {
                node_rX$initial_coop = node_rX$coop
                } else {
                node_rX$initial_coop = node_rX$initial_coop
                }
              node_rX$cost = (-50)*node_rX$coop*node_rX$prev_degree
              node_rX$n_coop_received = NA
              for (i in 1:(dim(node_rX)[1]))
                {
                node_rX[i,]$n_coop_received = sum(node_rX[node_rX$ego_id %in% node_rX[i,colnames(node_rX)[substr(colnames(node_rX),1,6) ==
                "alt_id"]],"coop"])
                }
              node_rX$benefit = 100*node_rX$n_coop_received
              node_rX$payoff = node_rX$cost + node_rX$benefit
              node_rX$wealth = node_rX$prev_wealth + node_rX$payoff
              node_rX$rel_rank = NA
              node_rX$local_rate_coop = NA
              for (i in 1:dim(node_rX)[1])
                {
                node_rX[i,]$rel_rank = rank1(na.omit(node_rX[node_rX$ego_id %in% node_rX[i,colnames(node_rX)[substr(colnames(node_rX),1,6) %in%
                c("ego_id","alt_id")]],"wealth"]))/length(na.omit(node_rX[node_rX$ego_id %in%
                node_rX[i,colnames(node_rX)[substr(colnames(node_rX),1,6) %in% c("ego_id","alt_id")]],"wealth"]))
                node_rX[i,]$local_rate_coop = mean(na.omit(node_rX[node_rX$ego_id %in% node_rX[i,colnames(node_rX)[substr(colnames(node_rX),1,6) %in%
                c("ego_id","alt_id")]],"coop"]))
                }
              node_rX$growth = as.numeric((node_rX$wealth/node_rX$prev_wealth) > 1)
              node_rX = node_rX[,c("ego_id","round","group","prev_degree","initial_wealth","initial_local_gini","initial_coop","coop","wealth","rel_rank","local_rate_coop","growth","everIsolated","maxDegreeLost")] #Pruning the previous-round data (degree is not updating yet)
              
              #3-2: Rewiring phase
              # 30% of ties (unidirectional) are being rewired
              link_rX_1 = link_import #Importing data (bidirectioanl ego-alter [ego_id < alter_id])
              colnames(link_rX_1) = c("ego_id","alt_id","prev_connected")
              link_rX_1$challenge = sample(0:1,dim(link_rX_1)[1],replace=TRUE,prob=c(0.7,0.3)) # The bidirectional ties being rewired are selected (rewiring rate = 0.3).
              ego_node_data =
              node_rX[,c("ego_id","wealth","coop","prev_degree","initial_wealth","initial_local_gini","initial_coop","rel_rank","local_rate_coop","growth")]
              colnames(ego_node_data) =
              c("ego_id","ego_wealth","ego_coop","ego_prev_degree","ego_initial_wealth","ego_initial_local_gini","ego_initial_coop","ego_rel_rank","ego_local_rate_coop","ego_growth")
              alt_node_data =
              node_rX[,c("ego_id","wealth","coop","prev_degree","initial_wealth","initial_local_gini","initial_coop","rel_rank","local_rate_coop","growth")]
              colnames(alt_node_data) =
              c("alt_id","alt_wealth","alt_coop","alt_prev_degree","alt_initial_wealth","alt_initial_local_gini","alt_initial_coop","alt_rel_rank","alt_local_rate_coop","alt_growth")
              link_rX_2 = merge(x=link_rX_1,y=ego_node_data,all.x=TRUE,all.y=FALSE,by="ego_id")
              link_rX_3 = merge(x=link_rX_2,y=alt_node_data,all.x=TRUE,all.y=FALSE,by="alt_id")
              link_rX_3$choice = sample(c("ego","alt"),dim(link_rX_3)[1],replace=TRUE,prob=c(0.5,0.5)) #decision maker for breaking a link, which is a unilateral decision
              #ego_prob: probability of choosing to connect when challenged (asked)
              
              if(modelForPrediction=="linear"){
                link_rX_3$ego_prob = inv.logit((0.5134401)*link_rX_3$prev_connected + (-0.852406)*link_rX_3$ego_coop + (2.96549)*link_rX_3$alt_coop + (-0.1808545))
                link_rX_3$alt_prob = inv.logit((0.5134401)*link_rX_3$prev_connected + (-0.852406)*link_rX_3$alt_coop + (2.96549)*link_rX_3$ego_coop + (-0.1808545))}
              if(modelForPrediction=="random forest"){
                link_rX_3$ego_prob = predict(model3,
                                             newdata=
                                               data.frame(
                                                 previouslyconnected = link_rX_3$prev_connected,
                                                 ego_behavior = link_rX_3$ego_coop,
                                                 alter_behavior = link_rX_3$alt_coop
                                                 ),
                                            type = "prob"
                                             )[[1]]$C
                link_rX_3$alt_prob = predict(model3,
                                             newdata=
                                               data.frame(
                                                 previouslyconnected = link_rX_3$prev_connected,
                                                 ego_behavior = link_rX_3$alt_coop,
                                                 alter_behavior = link_rX_3$ego_coop
                                                 ),
                                            type = "prob"
                                             )[[1]]$C
              }
              link_rX_3$prob_connect = ifelse(link_rX_3$prev_connected == 1, ifelse(link_rX_3$choice == "ego", link_rX_3$ego_prob,
              link_rX_3$alt_prob), link_rX_3$ego_prob*link_rX_3$alt_prob)
              link_rX_3$connect_update = apply(data.frame(link_rX_3$prob_connect),1, function(x) {sample(1:0,1,prob=c(x,(1-x)))})
              link_rX_3$connected = ifelse(link_rX_3$challenge==0,link_rX_3$prev_connected,link_rX_3$connect_update)
              link_rX = link_rX_3[,c("ego_id","alt_id","connected")] #pruning and data is updated
              #Reflect the degree and local gini coefficient into the node data
              link_rXc_ego = link_rX[link_rX$connected==1,]
              link_rXc_alt = link_rX[link_rX$connected==1,]
              colnames(link_rXc_alt) = c("alt_id","ego_id","connected")
              link_rXc = rbind(link_rXc_ego,link_rXc_alt)
              link_rXc = link_rXc[order(link_rXc$ego_id),]
              link_rXc$alternumber = NA
              link_rXc[1,]$alternumber = 1
              for (i in 1:(dim(link_rXc)[1]-1))
                {
                  if (link_rXc[i,]$ego_id == link_rXc[i+1,]$ego_id)
                  {
                  link_rXc[i+1,]$alternumber = link_rXc[i,]$alternumber + 1
                  }
                  else
                  {
                  link_rXc[i+1,]$alternumber = 1
                  }
                  #print(i)
                }
              link_rXc2 = reshape(link_rXc, direction = "wide", idvar=c("ego_id","connected"), timevar="alternumber")
              if(class(link_rXc2[,colnames(link_rXc2)[substr(colnames(link_rXc2),1,3) == "alt"]]) == "data.frame") {
                link_rXc2$degree = apply(link_rXc2[,colnames(link_rXc2)[substr(colnames(link_rXc2),1,3) == "alt"]],1,function(x) {length(na.omit(x))})
                }
              if(class(link_rXc2[,colnames(link_rXc2)[substr(colnames(link_rXc2),1,3) == "alt"]]) == "integer") {
                link_rXc2$degree = sapply(link_rXc2[,colnames(link_rXc2)[substr(colnames(link_rXc2),1,3) == "alt"]],function(x) {length(na.omit(x))})
                }
              node_rX_final = merge(x=node_rX[,c("ego_id","round","group","initial_wealth","initial_local_gini","initial_coop","coop","wealth","growth","everIsolated","maxDegreeLost")],y=link_rXc2,all.x=TRUE,all.y=FALSE,by="ego_id")
              node_rX_final[is.na(node_rX_final$degree)==1,"degree"] = 0
              node_rX_final$avg_env_wealth = NA
              node_rX_final$local_gini = NA #needs to be updated because the social network changes at the rewiring phase
              node_rX_final$local_rate_coop = NA
              node_rX_final$rel_rank = NA
              for (i in 1:dim(node_rX_final)[1])
                {
                node_rX_final[i,]$avg_env_wealth = mean(na.omit(node_rX_final[node_rX_final$ego_id %in%
                node_rX_final[i,colnames(node_rX_final)[substr(colnames(node_rX_final),1,6) %in% c("ego_id","alt_id")]],"wealth"]))
                node_rX_final[i,]$local_gini = gini(na.omit(node_rX_final[node_rX_final$ego_id %in%
                node_rX_final[i,colnames(node_rX_final)[substr(colnames(node_rX_final),1,6) %in% c("ego_id","alt_id")]],"wealth"]))
                node_rX_final[i,]$local_rate_coop = mean(na.omit(node_rX_final[node_rX_final$ego_id %in%
                node_rX_final[i,colnames(node_rX_final)[substr(colnames(node_rX_final),1,6) %in% c("ego_id","alt_id")]],"coop"]))
                node_rX_final[i,]$rel_rank = rank1(na.omit(node_rX_final[node_rX_final$ego_id %in%
                node_rX_final[i,colnames(node_rX_final)[substr(colnames(node_rX_final),1,6) %in%
                c("ego_id","alt_id")]],"wealth"]))/length(na.omit(node_rX_final[node_rX_final$ego_id %in%
                node_rX_final[i,colnames(node_rX_final)[substr(colnames(node_rX_final),1,6) %in% c("ego_id","alt_id")]],"wealth"]))
                node_rX_final[i,]$everIsolated = ifelse(node_rX_final[i,]$everIsolated==1,1,ifelse(node_rX_final[i,]$degree<=isolationDegree,1,0))
                node_rX_final[i,]$maxDegreeLost = pmax(node_r0[i,]$initial_degree - node_rX_final[i,]$degree, node_rX_final[i,]$maxDegreeLost, na.rm=TRUE)
              }
              
              
              #Finalization of round X and Visualization
              #plot(graph.data.frame(link_rX[link_rX$connected==1,],directed=F)) #plot.igraph
              result[result$round==k,2:25] =
              c(length(node_rX_final$ego_id),length(node_rX_final[node_rX_final$group=="rich",]$ego_id),mean(node_rX_final$coop),mean(node_rX_final$degree),mean(node_rX_final$wealth),gini(node_rX_final$wealth),gmd(node_rX_final$wealth),mean(node_rX_final[node_rX_final$group=="rich",]$coop),mean(node_rX_final[node_rX_final$group=="rich",]$degree),mean(node_rX_final[node_rX_final$group=="rich",]$wealth),gini(node_rX_final[node_rX_final$group=="rich",]$wealth),gmd(node_rX_final[node_rX_final$group=="rich",]$wealth),mean(node_rX_final[node_rX_final$group=="poor",]$coop),mean(node_rX_final[node_rX_final$group=="poor",]$degree),mean(node_rX_final[node_rX_final$group=="poor",]$wealth),gini(node_rX_final[node_rX_final$group=="poor",]$wealth),gmd(node_rX_final[node_rX_final$group=="poor",]$wealth),
                                             as.numeric(ifelse(is.na(table(node_rX_final$degree<=isolationDegree)["TRUE"]),0,1)),
                                             as.numeric(sum(node_rX_final$everIsolated)/length(node_rX_final$ego_id)),
                prop.table(table(node_rX_final[node_rX_final$initial_coop==1]$coop))["0"],
                prop.table(table(node_rX_final[node_rX_final$initial_coop==0]$coop))["1"],
                suppressWarnings({mean(node_rX_final$maxDegreeLost,na.rm=TRUE)}),
                suppressWarnings({mean(node_rX_final[node_rX_final$initial_coop==1]$maxDegreeLost,na.rm=TRUE)}),
                suppressWarnings({mean(node_rX_final[node_rX_final$initial_coop==0]$maxDegreeLost,na.rm=TRUE)})
                )
              
              #For the loop
              node_import = node_rX_final
              colnames(node_import)[colnames(node_import) %in%
              c("coop","wealth","growth","degree","avg_env_wealth","local_gini","local_rate_coop","rel_rank")] =
              c("prev_coop","prev_wealth","prev_growth","prev_degree","prev_avg_env_wealth","prev_local_gini","prev_local_rate_coop","prev_rel_rank")
              link_import = link_rX
              #print(paste0("Round ",k," is done."))
            }
            
           
            trends.df = rbind(trends.df,cbind(result[c("round","gini","gmd","avg_wealth","avg_coop","avg_degree")],V,GINI,fractionCoop))
            
            link_rX_final = data.table::melt(setDT(node_rX_final),
                                     measure = patterns('alt_id'),
                                     variable.name = 'linkNumber', 
                                     value.name = c('alt_id'))
            link_rX_final = data.frame(link_rX_final)[c("ego_id","alt_id")]
            link_rX_final = link_rX_final[complete.cases(link_rX_final),]
            link_rX_final = data.frame(t(unique(apply(link_rX_final, 1, function(x) sort(x))))) %>% distinct(X1, X2)
            node_g_final = data.frame(node_rX_final)[c("ego_id","initial_coop","coop")]
            node_g_final$initial_coop = factor(node_g_final$initial_coop)
            
            g_rX_final = graph_from_data_frame(link_rX_final, directed = FALSE, vertices=node_g_final)
            g_r0 = graph_from_data_frame(link_r0[link_r0$connected==1,][1:2], directed = FALSE, vertices=node_r0)
            
            E(g_r0)$coopEdgeC = sapply(E(g_r0), function(e) prod(ifelse(V(g_r0)[inc(e)]$coop_rp=="C",1,0)))
            E(g_r0)$coopEdgeD = sapply(E(g_r0), function(e) prod(ifelse(V(g_r0)[inc(e)]$coop_rp=="D",1,0)))
            E(g_r0)$coopEdgeCD = sapply(E(g_r0), function(e) ifelse(sum(ifelse(V(g_r0)[inc(e)]$coop_rp=="C",1,0))==1,1,0))
            
            #C-assortativity, defined as number of observed C-C edges out of total possible C-C edges
            homophilyC[m] = sum(E(g_r0)$coopEdgeC) / (table(V(g_r0)$coop_rp)["C"]*(table(V(g_r0)$coop_rp)["C"]-1)/2)
            #D-assortativity, defined as number of observed C-C edges out of total possible C-C edges
            homophilyD[m] = sum(E(g_r0)$coopEdgeD) / (table(V(g_r0)$coop_rp)["D"]*(table(V(g_r0)$coop_rp)["D"]-1)/2)
            #heterophily, defined as number of observed C-D edges out of total possible C-D edges
            heterophily[m] = sum(E(g_r0)$coopEdgeCD) / (table(V(g_r0)$coop_rp)["C"]*table(V(g_r0)$coop_rp)["D"])
            
            coopFrac[m] = fractionCoop
            avgCoop[m] = prop.table(table(V(g_r0)$coop_rp))["C"]
            avgCoopFinal[m] = result[result$round==10,]$avg_coop
            percentIsolation[m] = max(result[result$round>=1,]$percentIsolation)
            isolationPersonRounds[m] = sum(N*result[result$round>=1,]$percentIsolation)
            isolation[m] = max(result[result$round>=1,]$isolation)
            #percentage of isolation among those who cooperated in both practice rounds
            percentIsolationC[m] = sum(node_rX_final[coop_rp_init==1,]$everIsolated)/length(node_rX_final[coop_rp_init==1,]$everIsolated)
            #percentage of isolation among those who defected at least once in practice rounds
            percentIsolationD[m] = sum(node_rX_final[coop_rp_init<=0.5,]$everIsolated)/length(node_rX_final[coop_rp_init<=0.5,]$everIsolated)
            
            nCommunities[m] = max(membership(cluster_louvain(g_rX_final)),na.rm=TRUE)
            communitySize[m] = mean(table(membership(cluster_louvain(g_rX_final))),na.rm=TRUE)
            assortativityInitial[m] = assortativity(g_r0, V(g_r0)$coop_rp == "C")
            assortativityFinal[m] = assortativity(g_rX_final, V(g_r0)$coop_rp == "C")
            conversionRate[m] = prop.table(table(V(g_rX_final)$coop == ifelse(V(g_r0)$coop_rp=="C","1","0")))["FALSE"]
            conversionToD[m] = prop.table(table(V(g_rX_final)$coop[V(g_r0)$coop_rp == "C"]))["0"]
            conversionToC[m] = prop.table(table(V(g_rX_final)$coop[V(g_r0)$coop_rp == "C"]))["1"]
            transitivity[m] = mean(transitivity(g_rX_final, type="global"),na.rm=TRUE)
            degree[m] = mean(igraph::degree(g_rX_final),na.rm=TRUE)
            degreeC[m] = mean(igraph::degree(g_r0)[coop_rp_init==1],na.rm=TRUE)
            degreeD[m] = mean(igraph::degree(g_r0)[coop_rp_init<=0.5],na.rm=TRUE)
            meanConversionToD[m] = mean(result[result$round>=2,]$meanConversionToD, na.rm=TRUE)
            meanConversionToC[m] = mean(result[result$round>=2,]$meanConversionToC, na.rm=TRUE)
            degreeLost[m] = result[result$round==10,]$degreeLost
            degreeLostC[m] = result[result$round==10,]$degreeLostC
            degreeLostD[m] = result[result$round==10,]$degreeLostD
            avg_wealth[m] = result[result$round==10,]$avg_wealth
            gini[m] = result[result$round==10,]$gini
        
          }
          
        df.netIntLowDegree = rbind(df.netIntLowDegree,
                          data.frame(
                            coopFrac = coopFrac,
                            avgCoop = avgCoop,
                            avgCoopFinal = avgCoopFinal,
                            percentIsolation = percentIsolation,
                            isolationPersonRounds = isolationPersonRounds,
                            isolation = isolation,
                            percentIsolationC = percentIsolationC,
                            percentIsolationD = percentIsolationD,
                            nCommunities = nCommunities,
                            communitySize = communitySize,
                            assortativityInitial = assortativityInitial,
                            assortativityFinal = assortativityFinal,
                            conversionRate = conversionRate,
                            conversionToD = conversionToD, 
                            conversionToC = conversionToC, 
                            homophilyC = homophilyC,
                            homophilyD = homophilyD,
                            heterophily = heterophily,
                            transitivity = transitivity, 
                            degree = degree, 
                            degreeC = degreeC, 
                            degreeD = degreeD,
                            meanConversionToD = meanConversionToD, 
                            meanConversionToC = meanConversionToC,
                            degreeLost = degreeLost,
                            degreeLostC = degreeLostC,
                            degreeLostD = degreeLostD,
                            avg_wealth = avg_wealth,
                            gini = gini
                            ))
        #plot(g_r0,vertex.color=V(g_rX_final)$initial_coop,vertex.label=ifelse(is.na(V(g_rX_final)$initial_coop),"NA",ifelse(V(g_rX_final)$initial_coop==1,"C","D")),main=paste("fracCoop=",frac,", round 0",sep=""))
        #plot(g_rX_final,vertex.color=V(g_rX_final)$initial_coop,vertex.label=ifelse(is.na(V(g_rX_final)$initial_coop),"NA",ifelse(V(g_rX_final)$initial_coop==1,"C","D")),main=paste("fracCoop=",frac,", final round",sep=""))
        
      }
      
    sum.netIntLowDegree <- data.frame(
      df.netIntLowDegree %>% 
        group_by(coopFrac) %>%
          summarise(
            mean.isolation = mean(isolation),
            ci.isolation   = 1.96 * sd(isolation)/sqrt(n()),
            mean.percentIsolation = mean(percentIsolation),
            ci.percentIsolation   = 1.96 * sd(percentIsolation)/sqrt(n()),
            mean.isolationPersonRounds = mean(isolationPersonRounds),
            ci.isolationPersonRounds   = 1.96 * sd(isolationPersonRounds)/sqrt(n()),
            mean.percentIsolationC = mean(percentIsolationC,na.rm=TRUE),
            ci.percentIsolationC   = 1.96 * sd(percentIsolationC,na.rm=TRUE)/sqrt(sum(isolation)),
            mean.percentIsolationD = mean(percentIsolationD,na.rm=TRUE),
            ci.percentIsolationD   = 1.96 * sd(percentIsolationD,na.rm=TRUE)/sqrt(sum(isolation)),
            mean.avgCoop = mean(avgCoop,na.rm=TRUE),
            ci.avgCoop   = 1.96 * sd(avgCoop,na.rm=TRUE)/sqrt(n()),
            mean.avgCoopFinal = mean(avgCoopFinal,na.rm=TRUE),
            ci.avgCoopFinal   = 1.96 * sd(avgCoopFinal,na.rm=TRUE)/sqrt(n()),
            mean.nCommunities = mean(nCommunities,na.rm=TRUE),
            ci.nCommunities   = 1.96 * sd(nCommunities,na.rm=TRUE)/sqrt(n()),
            mean.communitySize = mean(communitySize,na.rm=TRUE),
            ci.communitySize   = 1.96 * sd(communitySize,na.rm=TRUE)/sqrt(n()),
            mean.assortativityInitial = mean(assortativityInitial,na.rm=TRUE),
            ci.assortativityInitial   = 1.96 * sd(assortativityInitial,na.rm=TRUE)/sqrt(n()),
            mean.assortativityFinal = mean(assortativityFinal,na.rm=TRUE),
            ci.assortativityFinal   = 1.96 * sd(assortativityFinal,na.rm=TRUE)/sqrt(n()),
            mean.conversionRate = mean(conversionRate,na.rm=TRUE),
            ci.conversionRate   = 1.96 * sd(conversionRate,na.rm=TRUE)/sqrt(n()),
            mean.conversionToD = mean(conversionToD,na.rm=TRUE),
            ci.conversionToD   = 1.96 * sd(conversionToD,na.rm=TRUE)/sqrt(n()),
            mean.conversionToC = mean(conversionToC,na.rm=TRUE),
            ci.conversionToC   = 1.96 * sd(conversionToC,na.rm=TRUE)/sqrt(n()),
            mean.homophilyC = mean(homophilyC,na.rm=TRUE),
            ci.homophilyC   = 1.96 * sd(homophilyC,na.rm=TRUE)/sqrt(n()),
            mean.homophilyD = mean(homophilyD,na.rm=TRUE),
            ci.homophilyD   = 1.96 * sd(homophilyD,na.rm=TRUE)/sqrt(n()),
            mean.heterophily = mean(heterophily,na.rm=TRUE),
            ci.heterophily   = 1.96 * sd(heterophily,na.rm=TRUE)/sqrt(n()),
            mean.transitivity = mean(transitivity,na.rm=TRUE),
            ci.transitivity   = 1.96 * sd(transitivity,na.rm=TRUE)/sqrt(n()),
            mean.degree = mean(degree,na.rm=TRUE),
            ci.degree   = 1.96 * sd(degree,na.rm=TRUE)/sqrt(n()),
            mean.degreeC = mean(degreeC,na.rm=TRUE),
            ci.degreeC   = 1.96 * sd(degreeC,na.rm=TRUE)/sqrt(n()),
            mean.degreeD = mean(degreeD,na.rm=TRUE),
            ci.degreeD   = 1.96 * sd(degreeD,na.rm=TRUE)/sqrt(n()),
            mean.meanConversionToD = mean(meanConversionToD,na.rm=TRUE),
            ci.meanConversionToD   = 1.96 * sd(meanConversionToD,na.rm=TRUE)/sqrt(n()),
            mean.meanConversionToC = mean(meanConversionToC,na.rm=TRUE),
            ci.meanConversionToC   = 1.96 * sd(meanConversionToC,na.rm=TRUE)/sqrt(n()),
            mean.degreeLost = mean(degreeLost,na.rm=TRUE),
            ci.degreeLost   = 1.96 * sd(degreeLost,na.rm=TRUE)/sqrt(n()),
            mean.degreeLostC = mean(degreeLostC,na.rm=TRUE),
            ci.degreeLostC   = 1.96 * sd(degreeLostC,na.rm=TRUE)/sqrt(n()),
            mean.degreeLostD = mean(degreeLostD,na.rm=TRUE),
            ci.degreeLostD   = 1.96 * sd(degreeLostD,na.rm=TRUE)/sqrt(n()),
            mean.avg_wealth = mean(avg_wealth,na.rm=TRUE),
            ci.avg_wealth   = 1.96 * sd(avg_wealth,na.rm=TRUE)/sqrt(n()),
            mean.gini = mean(gini,na.rm=TRUE),
            ci.gini   = 1.96 * sd(gini,na.rm=TRUE)/sqrt(n())
            )
      )
          
    kable(sum.netIntLowDegree[c(1:9)]) %>% kableExtra::kable_styling(font_size = 10)
    kable(sum.netIntLowDegree[c(1,10:17)]) %>% kableExtra::kable_styling(font_size = 10) 
    kable(sum.netIntLowDegree[c(1,18:25)]) %>% kableExtra::kable_styling(font_size = 10) 
    kable(sum.netIntLowDegree[c(1,26:33)]) %>% kableExtra::kable_styling(font_size = 10) 
    kable(sum.netIntLowDegree[c(1,34:ncol(sum.netIntLowDegree))]) %>% kableExtra::kable_styling(font_size = 10) 
    
    compare_means(percentIsolation ~ coopFrac, data=df.netIntLowDegree)
    compare_means(isolationPersonRounds ~ coopFrac, data=df.netIntLowDegree)
    compare_means(avgCoop ~ coopFrac, data=df.netIntLowDegree)
    compare_means(avgCoopFinal ~ coopFrac, data=df.netIntLowDegree)
    compare_means(nCommunities ~ coopFrac, data=df.netIntLowDegree)
    compare_means(communitySize ~ coopFrac, data=df.netIntLowDegree)
    compare_means(assortativityInitial ~ coopFrac, data=df.netIntLowDegree)
    compare_means(assortativityFinal ~ coopFrac, data=df.netIntLowDegree)
    #compare_means(conversionRate ~ coopFrac, data=df.netIntLowDegree)
    #compare_means(conversionToD ~ coopFrac, data=df.netIntLowDegree)
    #compare_means(conversionToC ~ coopFrac, data=df.netIntLowDegree)
    #compare_means(degreeC ~ coopFrac, data=df.netIntLowDegree)
    #compare_means(degreeD ~ coopFrac, data=df.netIntLowDegree)
    #compare_means(meanConversionToD ~ coopFrac, data=df.netIntLowDegree)
    #compare_means(meanConversionToC ~ coopFrac, data=df.netIntLowDegree)
    #compare_means(degreeLost ~ coopFrac, data=df.netIntLowDegree)
    #compare_means(degreeLostC ~ coopFrac, data=df.netIntLowDegree)
    #compare_means(degreeLostD ~ coopFrac, data=df.netIntLowDegree)
    
    summary(lm(percentIsolation ~ assortativityInitial, data=df.netIntLowDegree))
    
    #plot(df.netIntLowDegree$assortativityInitial, df.netIntLowDegree$percentIsolation)
    
    
    
    #percentIsolation
    g.percentIsolation = ggbarplot(data=df.netIntLowDegree, x="coopFrac", y="percentIsolation", add = "mean_se", color="coopFrac") +
      stat_compare_means(ref.group = "0", label = "p.signif", label.y = 0.098, method="t.test", color="black") +  
      labs(
        title = paste("Isolation when defectors are assigned to 25% of nodes by degree, ","V=",V,", Gini=",GINI,sep=""),
        x = "Degree percentile of nodes assigned to defectors ",
        y = "Propoption of ever-isolated individuals") +
      annotate("text", x=1, y=0.0990, label= "ref", color="black") +
      annotate("text", x=2.4, y= -0.0022, label= "Lowest degree nodes assigned to defectors", size=2.5) +
      geom_segment(aes(x = 3.3, y = -0.0024, xend = 3.7, yend = -0.0024), linewidth=0.2, arrow = arrow(length = unit(0.1, "cm"))) +
      annotate("text", x=4.6, y= -0.0022, label= "Highest degree nodes assigned to defectors", size=2.5) +
      theme_bw() +
      theme(plot.title = element_text(hjust = 0.5, size=12),legend.position="none") +
      coord_cartesian(ylim=c(0,0.10)) + 
      scale_x_discrete(labels=c('Control','0-25','25-50','50-75','75-100')) +
      scale_color_manual(values = c('0' = "black",'0.25'="black",'0.5'="black",'0.75'="black",'1'="black")) +
      geom_vline(xintercept = 1.5, linetype = "longdash")
    print(g.percentIsolation)
    
    #isolationPersonRounds
    g.isolationPersonRounds = ggbarplot(data=df.netIntLowDegree, x="coopFrac", y="isolationPersonRounds", add = "mean_se", color="coopFrac") +
      stat_compare_means(ref.group = "0", label = "p.signif", label.y = 5.98, method="t.test", color="black") +  
      labs(
        title = paste("Isolation when defectors are assigned to 25% of nodes by degree, ","V=",V,", Gini=",GINI,sep=""),
        x = "Degree percentile of nodes assigned to defectors ",
        y = "Person-rounds of isolation") +
      annotate("text", x=1, y=6.080, label= "ref", color="black") +
      annotate("text", x=2.4, y= -0.15, label= "Lowest degree nodes assigned to defectors", size=2.5) +
      geom_segment(aes(x = 3.3, y = -0.17, xend = 3.7, yend = -0.17), linewidth=0.2, arrow = arrow(length = unit(0.1, "cm"))) +
      annotate("text", x=4.6, y= -0.15, label= "Highest degree nodes assigned to defectors", size=2.5) +
      theme_bw() +
      theme(plot.title = element_text(hjust = 0.5, size=12),legend.position="none") +
      coord_cartesian(ylim=c(0,6)) + 
      scale_x_discrete(labels=c('Control','0-25','25-50','50-75','75-100')) +
      scale_color_manual(values = c('0' = "black",'0.25'="black",'0.5'="black",'0.75'="black",'1'="black")) +
      geom_vline(xintercept = 1.5, linetype = "longdash")
    print(g.isolationPersonRounds)
    
    #percentIsolationC
    #percentage of isolation among those who cooperated in both practice rounds
    g.percentIsolationC = ggbarplot(data=df.netIntLowDegree, x="coopFrac", y="percentIsolationC", add = "mean_se") +
      stat_compare_means(ref.group = "0", label = "p.signif", label.y = 0.098, method="t.test", color="black") +  
      labs(
        title = paste("Isolation among initial cooperators, ","V=",V,", Gini=",GINI,sep=""),
        x = "Degree percentile of nodes assigned to defectors ",
        y = "Propoption of ever-isolated individuals") +
      annotate("text", x=1, y=0.0990, label= "ref", color="black") +
      annotate("text", x=2.4, y= -0.0022, label= "Lowest degree nodes assigned to defectors", size=2.5) +
      geom_segment(aes(x = 3.3, y = -0.0024, xend = 3.7, yend = -0.0024), linewidth=0.2, arrow = arrow(length = unit(0.1, "cm"))) +
      annotate("text", x=4.6, y= -0.0022, label= "Highest degree nodes assigned to defectors", size=2.5) +
      theme_bw() +
      theme(plot.title = element_text(hjust = 0.5, size=12),legend.position="none") +
      coord_cartesian(ylim=c(0,0.10)) + 
      scale_x_discrete(labels=c('Control','0-25','25-50','50-75','75-100')) +
      scale_color_manual(values = c('0' = "black",'0.25'="black",'0.5'="black",'0.75'="black",'1'="black")) +
      geom_vline(xintercept = 1.5, linetype = "longdash")
    print(g.percentIsolationC)
    
    #percentIsolationD
    #percentage of isolation among those who defected at least once in practice rounds
    g.percentIsolationD = ggbarplot(data=df.netIntLowDegree, x="coopFrac", y="percentIsolationD", add = "mean_se") +
      stat_compare_means(ref.group = "0", label = "p.signif", label.y = 0.298, method="t.test", color="black") +  
      labs(
        title = paste("Isolation among initial defectors, ","V=",V,", Gini=",GINI,sep=""),
        x = "Degree percentile of nodes assigned to defectors ",
        y = "Propoption of ever-isolated individuals") +
      annotate("text", x=1, y=0.2990, label= "ref", color="black") +
      annotate("text", x=2.4, y= -0.0062, label= "Lowest degree nodes assigned to defectors", size=2.5) +
      geom_segment(aes(x = 3.3, y = -0.0064, xend = 3.7, yend = -0.0064), linewidth=0.2, arrow = arrow(length = unit(0.1, "cm"))) +
      annotate("text", x=4.6, y= -0.0062, label= "Highest degree nodes assigned to defectors", size=2.5) +
      theme_bw() +
      theme(plot.title = element_text(hjust = 0.5, size=12),legend.position="none") +
      coord_cartesian(ylim=c(0,0.30)) + 
      scale_x_discrete(labels=c('Control','0-25','25-50','50-75','75-100')) +
      scale_color_manual(values = c('0' = "black",'0.25'="black",'0.5'="black",'0.75'="black",'1'="black")) +
      geom_vline(xintercept = 1.5, linetype = "longdash")
    print(g.percentIsolationD)
      
    #avgCoopFinal
    g.avgCoopFinal = ggbarplot(data=df.netIntLowDegree, x="coopFrac", y="avgCoopFinal", add = "mean_se") +
      stat_compare_means(ref.group = "0", label = "p.signif", label.y = 0.98, method="t.test", color="black") +  
      labs(
        title = paste("Cooperation in final round, ","V=",V,", Gini=",GINI,sep=""),
        x = "Degree percentile of nodes assigned to defectors ",
        y = "Propoption of cooperators in final round") +
      annotate("text", x=1, y=0.990, label= "ref", color="black") +
      annotate("text", x=2.4, y= -0.0212, label= "Lowest degree nodes assigned to defectors", size=2.5) +
      geom_segment(aes(x = 3.3, y = -0.0214, xend = 3.7, yend = -0.0214), linewidth=0.2, arrow = arrow(length = unit(0.1, "cm"))) +
      annotate("text", x=4.6, y= -0.0212, label= "Highest degree nodes assigned to defectors", size=2.5) +
      theme_bw() +
      theme(plot.title = element_text(hjust = 0.5, size=12),legend.position="none") +
      coord_cartesian(ylim=c(0,1.0)) + 
      scale_x_discrete(labels=c('Control','0-25','25-50','50-75','75-100')) +
      scale_color_manual(values = c('0' = "black",'0.25'="black",'0.5'="black",'0.75'="black",'1'="black")) +
      geom_vline(xintercept = 1.5, linetype = "longdash")
    print(g.avgCoopFinal)
    
    #avg_wealth
    g.avg_wealth = ggbarplot(data=df.netIntLowDegree, x="coopFrac", y="avg_wealth", add = "mean_se") +
      stat_compare_means(ref.group = "0", label = "p.signif", label.y = 6800, method="t.test", color="black") +  
      labs(
        title = paste("Wealth in final round, ","V=",V,", Gini=",GINI,sep=""),
        x = "Degree percentile of nodes assigned to defectors ",
        y = "Average wealth in final round") +
      annotate("text", x=1, y=6900, label= "ref", color="black") +
      annotate("text", x=2.4, y= -162, label= "Lowest degree nodes assigned to defectors", size=2.5) +
      geom_segment(aes(x = 3.3, y = -164, xend = 3.7, yend = -164), linewidth=0.2, arrow = arrow(length = unit(0.1, "cm"))) +
      annotate("text", x=4.6, y= -162, label= "Highest degree nodes assigned to defectors", size=2.5) +
      theme_bw() +
      theme(plot.title = element_text(hjust = 0.5, size=12),legend.position="none") +
      coord_cartesian(ylim=c(0,7000)) + 
      scale_x_discrete(labels=c('Control','0-25','25-50','50-75','75-100')) +
      scale_color_manual(values = c('0' = "black",'0.25'="black",'0.5'="black",'0.75'="black",'1'="black")) +
      geom_vline(xintercept = 1.5, linetype = "longdash")
    print(g.avg_wealth)
    
    #gini
    g.gini = ggbarplot(data=df.netIntLowDegree, x="coopFrac", y="gini", add = "mean_se") +
      stat_compare_means(ref.group = "0", label = "p.signif", label.y = 0.48, method="t.test", color="black") +  
      labs(
        title = paste("Gini coefficient in final round, ","V=",V,", Gini=",GINI,sep=""),
        x = "Degree percentile of nodes assigned to defectors ",
        y = "Gini coefficient in final round") +
      annotate("text", x=1, y=0.490, label= "ref", color="black") +
      annotate("text", x=2.4, y= -0.0112, label= "Lowest degree nodes assigned to defectors", size=2.5) +
      geom_segment(aes(x = 3.3, y = -0.0114, xend = 3.7, yend = -0.0114), linewidth=0.2, arrow = arrow(length = unit(0.1, "cm"))) +
      annotate("text", x=4.6, y= -0.0112, label= "Highest degree nodes assigned to defectors", size=2.5) +
      theme_bw() +
      theme(plot.title = element_text(hjust = 0.5, size=12),legend.position="none") +
      coord_cartesian(ylim=c(0,0.50)) + 
      scale_x_discrete(labels=c('Control','0-25','25-50','50-75','75-100')) +
      scale_color_manual(values = c('0' = "black",'0.25'="black",'0.5'="black",'0.75'="black",'1'="black")) +
      geom_vline(xintercept = 1.5, linetype = "longdash")
    print(g.gini)
    
    #degree
    g.degree = ggbarplot(data=df.netIntLowDegree, x="coopFrac", y="degree", add = "mean_se") +
      stat_compare_means(ref.group = "0", label = "p.signif", label.y = 14.8, method="t.test", color="black") +  
      labs(
        title = paste("Degree in final round, ","V=",V,", Gini=",GINI,sep=""),
        x = "Degree percentile of nodes assigned to defectors ",
        y = "Mean degree in final round") +
      annotate("text", x=1, y=14.90, label= "ref", color="black") +
      annotate("text", x=2.4, y= -0.312, label= "Lowest degree nodes assigned to defectors", size=2.5) +
      geom_segment(aes(x = 3.3, y = -0.314, xend = 3.7, yend = -0.314), linewidth=0.2, arrow = arrow(length = unit(0.1, "cm"))) +
      annotate("text", x=4.6, y= -0.312, label= "Highest degree nodes assigned to defectors", size=2.5) +
      theme_bw() +
      theme(plot.title = element_text(hjust = 0.5, size=12),legend.position="none") +
      coord_cartesian(ylim=c(0,15)) + 
      scale_x_discrete(labels=c('Control','0-25','25-50','50-75','75-100')) +
      scale_color_manual(values = c('0' = "black",'0.25'="black",'0.5'="black",'0.75'="black",'1'="black")) +
      geom_vline(xintercept = 1.5, linetype = "longdash")
    print(g.degree)
    
    #transitivity
    g.transitivity = ggbarplot(data=df.netIntLowDegree, x="coopFrac", y="transitivity", add = "mean_se") +
      stat_compare_means(ref.group = "0", label = "p.signif", label.y = 0.98, method="t.test", color="black") +  
      labs(
        title = paste("Transitivity in final round, ","V=",V,", Gini=",GINI,sep=""),
        x = "Degree percentile of nodes assigned to defectors ",
        y = "Transitivity in final round") +
      annotate("text", x=1, y=0.99, label= "ref", color="black") +
      annotate("text", x=2.4, y= -0.0212, label= "Lowest degree nodes assigned to defectors", size=2.5) +
      geom_segment(aes(x = 3.3, y = -0.0214, xend = 3.7, yend = -0.0214), linewidth=0.2, arrow = arrow(length = unit(0.1, "cm"))) +
      annotate("text", x=4.6, y= -0.0212, label= "Highest degree nodes assigned to defectors", size=2.5) +
      theme_bw() +
      theme(plot.title = element_text(hjust = 0.5, size=12),legend.position="none") +
      coord_cartesian(ylim=c(0,1.00)) + 
      scale_x_discrete(labels=c('Control','0-25','25-50','50-75','75-100')) +
      scale_color_manual(values = c('0' = "black",'0.25'="black",'0.5'="black",'0.75'="black",'1'="black")) +
      geom_vline(xintercept = 1.5, linetype = "longdash")
    print(g.transitivity)
    
    #initial C-assortativity     
    plotList <- lapply(
              unique(df.netIntLowDegree$coopFrac),
              function(key) {
                if(key==0){
                  ggplot(data = df.netIntLowDegree[df.netIntLowDegree$coopFrac==key,], aes(x = homophilyC, y = percentIsolation)) + 
                      geom_point() + 
                      scale_x_continuous(paste("C-assortativity, ","Control",sep="")) +
                      scale_y_continuous("Proportion isolated") +
                      geom_smooth(method='lm', formula= y~x) + 
                      stat_cor(method = "pearson")
                }
                else{
                  ggplot(data = df.netIntLowDegree[df.netIntLowDegree$coopFrac==key,], aes(x = homophilyC, y = percentIsolation)) + 
                      geom_point() + 
                      scale_x_continuous(paste("C-assortativity, degree %ile = ",key,sep="")) +
                      scale_y_continuous("Proportion isolated") +
                      geom_smooth(method='lm', formula= y~x) + 
                      stat_cor(method = "pearson")
                }
              }
    )
    plot= ggarrange(plotlist=plotList)
    print(annotate_figure(plot, top = text_grob(paste("Proportion of ever-isolated individuals, ","V=",V,", Gini=", GINI, sep=""), color = "black", face = "bold", size = 10)))
    
    lapply(unique(df.netIntLowDegree$coopFrac),
              function(key) {
                if(key==0){
                  reg = lm(percentIsolation ~ homophilyC + degreeD, data=df.netIntLowDegree[df.netIntLowDegree$coopFrac==key,])
                  print(paste("Regression on proportion of ever-isolated individuals, ","Control"," ; ",sep=""))
                  print(summary(reg)[4]$coefficients)
                }
                else{
                  reg = lm(percentIsolation ~ homophilyC + degreeD, data=df.netIntLowDegree[df.netIntLowDegree$coopFrac==key,])
                  print(paste("Regression on proportion of ever-isolated individuals, degree %ile = ",key," ; ",sep=""))
                  print(summary(reg)[4]$coefficients)
                }
              }
    )
    
    
    
    #initial D-assortativity     
    plotList <- lapply(
              unique(df.netIntLowDegree$coopFrac),
              function(key) {
                if(key==0){
                  ggplot(data = df.netIntLowDegree[df.netIntLowDegree$coopFrac==key,], aes(x = homophilyD, y = percentIsolation)) + 
                      geom_point() +
                      scale_x_continuous(paste("D-assortativity, ","Control",sep="")) + 
                      scale_y_continuous("Proportion isolated") +
                      geom_smooth(method='lm', formula= y~x) + 
                      stat_cor(method = "pearson")
                }
                else{
                  ggplot(data = df.netIntLowDegree[df.netIntLowDegree$coopFrac==key,], aes(x = homophilyD, y = percentIsolation)) + 
                      geom_point() +
                      scale_x_continuous(paste("D-assortativity, degree %ile = ",key,sep="")) + 
                      scale_y_continuous("Proportion isolated") +
                      geom_smooth(method='lm', formula= y~x) + 
                      stat_cor(method = "pearson")
                }
              }
    )
    plot= ggarrange(plotlist=plotList)
    print(annotate_figure(plot, top = text_grob(paste("Proportion of ever-isolated individuals, ","V=",V,", Gini=", GINI, sep=""), color = "black", face = "bold", size = 10)))
    
    lapply(unique(df.netIntLowDegree$coopFrac),
              function(key) {
                if(key==0){
                  reg = lm(percentIsolation ~ homophilyD + degreeD, data=df.netIntLowDegree[df.netIntLowDegree$coopFrac==key,])
                  print(paste("Regression on proportion of ever-isolated individuals, ","Control"," ; ",sep=""))
                  print(summary(reg)[4]$coefficients)
                }
                else{
                  reg = lm(percentIsolation ~ homophilyD + degreeD, data=df.netIntLowDegree[df.netIntLowDegree$coopFrac==key,])
                  print(paste("Regression on proportion of ever-isolated individuals, degree %ile = ",key," ; ",sep=""))
                  print(summary(reg)[4]$coefficients)
                }
              }
    )
    
    }
}

## Loading data last updated on  2023-01-20 20:37:15 
## Call model1.invisible(redo=TRUE) to update data.
## Loading data last updated on  2023-01-20 20:39:47 
## Call model2.invisible(redo=TRUE) to update data.
## Loading data last updated on  2023-01-20 22:02:59 
## Call model3(redo=TRUE) to update data.

## Warning: Removed 2 rows containing non-finite values (`stat_summary()`).

## Warning: Removed 2 rows containing non-finite values (`stat_compare_means()`).

## Warning: Removed 7 rows containing non-finite values (`stat_summary()`).

## Warning: Removed 7 rows containing non-finite values (`stat_compare_means()`).

## Warning: Removed 8 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 8 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 8 rows containing missing values (`geom_point()`).

## Warning: Removed 8 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 8 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 8 rows containing missing values (`geom_point()`).

## Warning: Removed 4 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 4 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 4 rows containing missing values (`geom_point()`).

## Warning: Removed 6 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 6 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 6 rows containing missing values (`geom_point()`).

## Warning: Removed 10 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 10 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 10 rows containing missing values (`geom_point()`).

## [1] "Regression on proportion of ever-isolated individuals, Control ; "
##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.19867820 0.015425964 12.879468 2.887548e-30
## homophilyC  -0.06415897 0.020975923 -3.058696 2.432895e-03
## degreeD     -0.05240828 0.006363278 -8.236051 6.360042e-15
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.25 ; "
##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.17029364 0.013074943 13.024427 8.768337e-31
## homophilyC  -0.03173874 0.025914348 -1.224756 2.216716e-01
## degreeD     -0.04743390 0.007738758 -6.129394 2.907025e-09
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.5 ; "
##                 Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.182345050 0.015463144 11.7922363 1.587187e-26
## homophilyC  -0.008446159 0.025584548 -0.3301274 7.415396e-01
## degreeD     -0.056180012 0.007752733 -7.2464784 3.808359e-12
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.75 ; "
##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.13077454 0.015991015  8.178001 9.211218e-15
## homophilyC   0.02830190 0.023533129  1.202641 2.300997e-01
## degreeD     -0.03728951 0.007117094 -5.239429 3.110379e-07
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 1 ; "
##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.15424385 0.018933674  8.146535 1.174418e-14
## homophilyC  -0.03381455 0.023664535 -1.428912 1.541207e-01
## degreeD     -0.02974242 0.007026136 -4.233112 3.108871e-05

## Warning: Removed 18 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 18 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 18 rows containing missing values (`geom_point()`).

## Warning: Removed 21 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 21 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 21 rows containing missing values (`geom_point()`).

## Warning: Removed 14 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 14 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 14 rows containing missing values (`geom_point()`).

## Warning: Removed 15 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 15 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 15 rows containing missing values (`geom_point()`).

## Warning: Removed 6 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 6 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 6 rows containing missing values (`geom_point()`).

## [1] "Regression on proportion of ever-isolated individuals, Control ; "
##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.19100505 0.016290556 11.724894 4.423334e-26
## homophilyD   0.01567285 0.024878442  0.629977 5.292252e-01
## degreeD     -0.06046064 0.009057203 -6.675420 1.327127e-10
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.25 ; "
##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.18476715 0.013829960 13.359919 9.752168e-32
## homophilyD  -0.03418777 0.023441788 -1.458411 1.458642e-01
## degreeD     -0.05416518 0.008372582 -6.469352 4.465804e-10
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.5 ; "
##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.19374985 0.015233690 12.718511 1.289978e-29
## homophilyD  -0.05714689 0.020900598 -2.734223 6.647098e-03
## degreeD     -0.05190122 0.008068478 -6.432591 5.332171e-10
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.75 ; "
##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.13644648 0.016480250  8.279393 5.016397e-15
## homophilyD  -0.02602873 0.020560318 -1.265969 2.065690e-01
## degreeD     -0.03206216 0.007467751 -4.293416 2.422417e-05
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 1 ; "
##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.15786434 0.018602376  8.486246 1.102142e-15
## homophilyD  -0.02536208 0.021670381 -1.170357 2.428150e-01
## degreeD     -0.03223958 0.007004176 -4.602908 6.230354e-06
## Loading data last updated on  2023-01-20 21:50:22 
## Call model1.visible(redo=TRUE) to update data.
## Loading data last updated on  2023-01-20 21:54:00 
## Call model2.visible(redo=TRUE) to update data.
## Loading data last updated on  2023-01-20 22:02:59 
## Call model3(redo=TRUE) to update data.

## Warning: Removed 2 rows containing non-finite values (`stat_summary()`).

## Warning: Removed 2 rows containing non-finite values (`stat_compare_means()`).

## Warning: Removed 10 rows containing non-finite values (`stat_summary()`).

## Warning: Removed 10 rows containing non-finite values (`stat_compare_means()`).

## Warning: Removed 5 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 5 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 5 rows containing missing values (`geom_point()`).

## Warning: Removed 6 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 6 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 6 rows containing missing values (`geom_point()`).

## Warning: Removed 9 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 9 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 9 rows containing missing values (`geom_point()`).

## Warning: Removed 5 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 5 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 5 rows containing missing values (`geom_point()`).

## Warning: Removed 8 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 8 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 8 rows containing missing values (`geom_point()`).

## [1] "Regression on proportion of ever-isolated individuals, Control ; "
##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.17571591 0.014860643 11.824246 1.359580e-26
## homophilyC  -0.06432891 0.022834910 -2.817130 5.178242e-03
## degreeD     -0.04505829 0.005827491 -7.732022 1.752765e-13
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.25 ; "
##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.13731226 0.012150360 11.301086 9.070973e-25
## homophilyC  -0.02947064 0.026301631 -1.120487 2.634364e-01
## degreeD     -0.03886852 0.006985811 -5.563924 6.015411e-08
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.5 ; "
##                Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.13751879 0.015005097  9.1648055 9.975156e-18
## homophilyC  -0.01400885 0.025078570 -0.5585985 5.768729e-01
## degreeD     -0.03923000 0.008126181 -4.8276064 2.253337e-06
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.75 ; "
##                Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.13716229 0.015989706  8.5781622 5.981214e-16
## homophilyC  -0.01985905 0.023429690 -0.8476021 3.973610e-01
## degreeD     -0.03486095 0.006827682 -5.1058251 5.985060e-07
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 1 ; "
##                Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.14184514 0.014152311 10.0227542 1.773029e-20
## homophilyC  -0.01757170 0.020125025 -0.8731268 3.833214e-01
## degreeD     -0.03218893 0.005194161 -6.1971377 1.984256e-09

## Warning: Removed 17 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 17 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 17 rows containing missing values (`geom_point()`).

## Warning: Removed 17 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 17 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 17 rows containing missing values (`geom_point()`).

## Warning: Removed 20 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 20 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 20 rows containing missing values (`geom_point()`).

## Warning: Removed 23 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 23 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 23 rows containing missing values (`geom_point()`).

## Warning: Removed 14 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 14 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 14 rows containing missing values (`geom_point()`).

## [1] "Regression on proportion of ever-isolated individuals, Control ; "
##                Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.15466757 0.014820199 10.4362679 9.393820e-22
## homophilyD   0.01159905 0.024748290  0.4686809 6.396624e-01
## degreeD     -0.04652123 0.008706247 -5.3434310 1.892087e-07
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.25 ; "
##                Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.14655674 0.012640012 11.5946679 1.194878e-25
## homophilyD  -0.02077082 0.021458917 -0.9679344 3.339125e-01
## degreeD     -0.04461620 0.007309784 -6.1036282 3.445650e-09
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.5 ; "
##                 Estimate  Std. Error     t value     Pr(>|t|)
## (Intercept)  0.145218425 0.015933440  9.11406582 1.652854e-17
## homophilyD  -0.001501206 0.020655615 -0.07267788 9.421149e-01
## degreeD     -0.044234605 0.007583763 -5.83280420 1.516757e-08
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.75 ; "
##                 Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.144989320 0.016705407  8.679185 3.622991e-16
## homophilyD   0.006690858 0.019505167  0.343030 7.318389e-01
## degreeD     -0.040511249 0.007255602 -5.583444 5.673009e-08
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 1 ; "
##                 Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.156797825 0.015268584 10.269310 3.092257e-21
## homophilyD   0.003402436 0.018823254  0.180757 8.566877e-01
## degreeD     -0.039357176 0.005804047 -6.780989 6.950120e-11
## Loading data last updated on  2023-01-20 20:37:15 
## Call model1.invisible(redo=TRUE) to update data.
## Loading data last updated on  2023-01-20 20:39:47 
## Call model2.invisible(redo=TRUE) to update data.
## Loading data last updated on  2023-01-20 22:02:59 
## Call model3(redo=TRUE) to update data.

## Warning: Removed 2 rows containing non-finite values (`stat_summary()`).

## Warning: Removed 2 rows containing non-finite values (`stat_compare_means()`).

## Warning: Removed 10 rows containing non-finite values (`stat_summary()`).

## Warning: Removed 10 rows containing non-finite values (`stat_compare_means()`).

## Warning: Removed 5 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 5 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 5 rows containing missing values (`geom_point()`).

## Warning: Removed 6 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 6 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 6 rows containing missing values (`geom_point()`).

## Warning: Removed 9 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 9 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 9 rows containing missing values (`geom_point()`).

## Warning: Removed 5 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 5 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 5 rows containing missing values (`geom_point()`).

## Warning: Removed 8 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 8 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 8 rows containing missing values (`geom_point()`).

## [1] "Regression on proportion of ever-isolated individuals, Control ; "
##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.18066201 0.016268071 11.105312 4.116412e-24
## homophilyC  -0.02661633 0.024997567 -1.064757 2.878717e-01
## degreeD     -0.04720626 0.006379403 -7.399792 1.481420e-12
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.25 ; "
##                 Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.173124014 0.012967160 13.3509580 5.347006e-32
## homophilyC   0.006150637 0.028069741  0.2191198 8.267114e-01
## degreeD     -0.054390584 0.007455427 -7.2954347 2.880397e-12
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.5 ; "
##                Estimate Std. Error   t value     Pr(>|t|)
## (Intercept)  0.17244380 0.01583961 10.886873 2.533305e-23
## homophilyC  -0.03479505 0.02647332 -1.314344 1.897835e-01
## degreeD     -0.04389057 0.00857812 -5.116572 5.714800e-07
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.75 ; "
##                Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.14190927 0.018725472  7.5784082 4.775540e-13
## homophilyC   0.01673700 0.027438403  0.6099847 5.423513e-01
## degreeD     -0.03345872 0.007995867 -4.1845020 3.795797e-05
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 1 ; "
##                Estimate Std. Error   t value     Pr(>|t|)
## (Intercept)  0.14605241 0.01631566  8.951672 4.422528e-17
## homophilyC  -0.03132559 0.02320137 -1.350161 1.780243e-01
## degreeD     -0.02786963 0.00598815 -4.654130 4.969476e-06

## Warning: Removed 17 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 17 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 17 rows containing missing values (`geom_point()`).

## Warning: Removed 17 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 17 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 17 rows containing missing values (`geom_point()`).

## Warning: Removed 20 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 20 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 20 rows containing missing values (`geom_point()`).

## Warning: Removed 23 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 23 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 23 rows containing missing values (`geom_point()`).

## Warning: Removed 14 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 14 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 14 rows containing missing values (`geom_point()`).

## [1] "Regression on proportion of ever-isolated individuals, Control ; "
##                Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.17839731 0.015874681 11.2378521 1.965734e-24
## homophilyD  -0.01234393 0.026509172 -0.4656473 6.418299e-01
## degreeD     -0.04876507 0.009325711 -5.2290991 3.338645e-07
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.25 ; "
##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.18745979 0.013379081 14.011410 3.669546e-34
## homophilyD  -0.03759836 0.022713633 -1.655321 9.897970e-02
## degreeD     -0.05450979 0.007737191 -7.045166 1.440067e-11
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.5 ; "
##                Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.17780231 0.016944240 10.4933776 6.572223e-22
## homophilyD  -0.01613386 0.021965984 -0.7344929 4.632697e-01
## degreeD     -0.04886891 0.008064868 -6.0594805 4.445981e-09
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.75 ; "
##                Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.15562793 0.019329491  8.0513207 2.516405e-14
## homophilyD  -0.01096610 0.022569037 -0.4858915 6.274325e-01
## degreeD     -0.03510733 0.008395312 -4.1817779 3.896512e-05
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 1 ; "
##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.15788640 0.016942337  9.319045 3.455588e-18
## homophilyD  -0.02948057 0.020886673 -1.411453 1.592094e-01
## degreeD     -0.03236170 0.006440291 -5.024882 8.941866e-07
## Loading data last updated on  2023-01-20 21:50:22 
## Call model1.visible(redo=TRUE) to update data.
## Loading data last updated on  2023-01-20 21:54:00 
## Call model2.visible(redo=TRUE) to update data.
## Loading data last updated on  2023-01-20 22:02:59 
## Call model3(redo=TRUE) to update data.

## Warning: Removed 2 rows containing non-finite values (`stat_summary()`).

## Warning: Removed 2 rows containing non-finite values (`stat_compare_means()`).

## Warning: Removed 10 rows containing non-finite values (`stat_summary()`).

## Warning: Removed 10 rows containing non-finite values (`stat_compare_means()`).

## Warning: Removed 5 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 5 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 5 rows containing missing values (`geom_point()`).

## Warning: Removed 6 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 6 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 6 rows containing missing values (`geom_point()`).

## Warning: Removed 9 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 9 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 9 rows containing missing values (`geom_point()`).

## Warning: Removed 5 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 5 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 5 rows containing missing values (`geom_point()`).

## Warning: Removed 8 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 8 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 8 rows containing missing values (`geom_point()`).

## [1] "Regression on proportion of ever-isolated individuals, Control ; "
##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.17070706 0.016570396 10.301930 2.059056e-21
## homophilyC  -0.05918856 0.025462121 -2.324573 2.078531e-02
## degreeD     -0.04107977 0.006497958 -6.321951 9.728131e-10
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.25 ; "
##                Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.15165120 0.012412076 12.2180360 5.871019e-28
## homophilyC  -0.01100552 0.026868162 -0.4096119 6.823941e-01
## degreeD     -0.04781037 0.007136284 -6.6996170 1.090793e-10
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.5 ; "
##                Estimate Std. Error   t value     Pr(>|t|)
## (Intercept)  0.13245604 0.01598670  8.285390 4.605159e-15
## homophilyC  -0.07375700 0.02671916 -2.760454 6.145107e-03
## degreeD     -0.02465288 0.00865778 -2.847483 4.725994e-03
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.75 ; "
##                Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.12903149 0.016365905  7.8841646 6.526496e-14
## homophilyC  -0.01359377 0.023980934 -0.5668573 5.712512e-01
## degreeD     -0.02949356 0.006988321 -4.2204069 3.268759e-05
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 1 ; "
##                Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.13763910 0.015138619  9.0919195 1.629815e-17
## homophilyC  -0.00496386 0.021527585 -0.2305814 8.178037e-01
## degreeD     -0.03026138 0.005556155 -5.4464607 1.102510e-07

## Warning: Removed 17 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 17 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 17 rows containing missing values (`geom_point()`).

## Warning: Removed 17 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 17 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 17 rows containing missing values (`geom_point()`).

## Warning: Removed 20 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 20 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 20 rows containing missing values (`geom_point()`).

## Warning: Removed 23 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 23 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 23 rows containing missing values (`geom_point()`).

## Warning: Removed 14 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 14 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 14 rows containing missing values (`geom_point()`).

## [1] "Regression on proportion of ever-isolated individuals, Control ; "
##                Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.14832377 0.016455568  9.0135918 3.211547e-17
## homophilyD  -0.02101279 0.027479197 -0.7646799 4.451061e-01
## degreeD     -0.03608365 0.009666958 -3.7326788 2.294423e-04
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.25 ; "
##                 Estimate Std. Error    t value     Pr(>|t|)
## (Intercept)  0.159161507 0.01279987 12.4346205 1.467049e-28
## homophilyD  -0.004760498 0.02173030 -0.2190719 8.267537e-01
## degreeD     -0.052903814 0.00740223 -7.1470108 7.714386e-12
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.5 ; "
##                Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.14433213 0.017151514  8.4151247 2.122536e-15
## homophilyD  -0.01839729 0.022234687 -0.8274141 4.087146e-01
## degreeD     -0.03851879 0.008163523 -4.7184026 3.775668e-06
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.75 ; "
##                 Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.137771466 0.016904719  8.149882 1.308941e-14
## homophilyD  -0.003850901 0.019737883 -0.195102 8.454577e-01
## degreeD     -0.033930159 0.007342169 -4.621272 5.878095e-06
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 1 ; "
##                 Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.151308846 0.015967152  9.476258 1.106920e-18
## homophilyD  -0.008860013 0.019684455 -0.450102 6.529815e-01
## degreeD     -0.034284891 0.006069593 -5.648631 3.937178e-08
## Loading data last updated on  2023-01-20 20:37:15 
## Call model1.invisible(redo=TRUE) to update data.
## Loading data last updated on  2023-01-20 20:39:47 
## Call model2.invisible(redo=TRUE) to update data.
## Loading data last updated on  2023-01-20 22:02:59 
## Call model3(redo=TRUE) to update data.

## Warning: Removed 2 rows containing non-finite values (`stat_summary()`).

## Warning: Removed 2 rows containing non-finite values (`stat_compare_means()`).

## Warning: Removed 10 rows containing non-finite values (`stat_summary()`).

## Warning: Removed 10 rows containing non-finite values (`stat_compare_means()`).

## Warning: Removed 5 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 5 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 5 rows containing missing values (`geom_point()`).

## Warning: Removed 6 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 6 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 6 rows containing missing values (`geom_point()`).

## Warning: Removed 9 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 9 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 9 rows containing missing values (`geom_point()`).

## Warning: Removed 5 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 5 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 5 rows containing missing values (`geom_point()`).

## Warning: Removed 8 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 8 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 8 rows containing missing values (`geom_point()`).

## [1] "Regression on proportion of ever-isolated individuals, Control ; "
##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.18850819 0.016884520 11.164557 2.583704e-24
## homophilyC  -0.04062833 0.025944805 -1.565953 1.184502e-01
## degreeD     -0.04751134 0.006621139 -7.175704 6.038767e-12
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.25 ; "
##                Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.18106291 0.013578601 13.3344306 6.133400e-32
## homophilyC  -0.02027101 0.029393312 -0.6896471 4.909695e-01
## degreeD     -0.05249748 0.007806973 -6.7244343 9.415237e-11
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.5 ; "
##                  Estimate  Std. Error     t value     Pr(>|t|)
## (Intercept)  0.1657495462 0.016513452 10.03724382 1.658447e-20
## homophilyC  -0.0009560503 0.027599540 -0.03464008 9.723909e-01
## degreeD     -0.0451592417 0.008943048 -5.04964742 7.895350e-07
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.75 ; "
##                Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.15409715 0.018310960  8.4155695 1.833500e-15
## homophilyC   0.01200799 0.026831019  0.4475415 6.548194e-01
## degreeD     -0.03668449 0.007818869 -4.6917898 4.184605e-06
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 1 ; "
##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.15625982 0.016538710  9.448127 1.245125e-18
## homophilyC  -0.04566776 0.023518557 -1.941776 5.313928e-02
## degreeD     -0.03056325 0.006070014 -5.035120 8.433722e-07

## Warning: Removed 17 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 17 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 17 rows containing missing values (`geom_point()`).

## Warning: Removed 17 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 17 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 17 rows containing missing values (`geom_point()`).

## Warning: Removed 20 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 20 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 20 rows containing missing values (`geom_point()`).

## Warning: Removed 23 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 23 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 23 rows containing missing values (`geom_point()`).

## Warning: Removed 14 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 14 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 14 rows containing missing values (`geom_point()`).

## [1] "Regression on proportion of ever-isolated individuals, Control ; "
##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.17813156 0.016500390 10.795597 6.045605e-23
## homophilyD  -0.03096313 0.027554046 -1.123723 2.620929e-01
## degreeD     -0.04444306 0.009693289 -4.584931 6.852697e-06
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.25 ; "
##                Estimate Std. Error    t value     Pr(>|t|)
## (Intercept)  0.18954846 0.01412758 13.4169131 4.938312e-32
## homophilyD  -0.02054521 0.02398435 -0.8566091 3.923935e-01
## degreeD     -0.05734412 0.00817005 -7.0188219 1.690763e-11
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.5 ; "
##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.17412849 0.017687624  9.844651 8.379347e-20
## homophilyD  -0.02560781 0.022929684 -1.116797 2.650488e-01
## degreeD     -0.04489274 0.008418693 -5.332507 2.013138e-07
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.75 ; "
##                 Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.159152504 0.018895764  8.4226551 2.094009e-15
## homophilyD  -0.002973771 0.022062619 -0.1347878 8.928786e-01
## degreeD     -0.036650867 0.008206933 -4.4658423 1.167045e-05
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 1 ; "
##                Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.17124919 0.017253333  9.9255716 4.064066e-20
## homophilyD  -0.01613857 0.021270072 -0.7587456 4.486362e-01
## degreeD     -0.03941305 0.006558509 -6.0094530 5.716298e-09
## Loading data last updated on  2023-01-20 21:50:22 
## Call model1.visible(redo=TRUE) to update data.
## Loading data last updated on  2023-01-20 21:54:00 
## Call model2.visible(redo=TRUE) to update data.
## Loading data last updated on  2023-01-20 22:02:59 
## Call model3(redo=TRUE) to update data.

## Warning: Removed 2 rows containing non-finite values (`stat_summary()`).

## Warning: Removed 2 rows containing non-finite values (`stat_compare_means()`).

## Warning: Removed 10 rows containing non-finite values (`stat_summary()`).

## Warning: Removed 10 rows containing non-finite values (`stat_compare_means()`).

## Warning: Removed 5 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 5 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 5 rows containing missing values (`geom_point()`).

## Warning: Removed 6 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 6 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 6 rows containing missing values (`geom_point()`).

## Warning: Removed 9 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 9 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 9 rows containing missing values (`geom_point()`).

## Warning: Removed 5 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 5 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 5 rows containing missing values (`geom_point()`).

## Warning: Removed 8 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 8 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 8 rows containing missing values (`geom_point()`).

## [1] "Regression on proportion of ever-isolated individuals, Control ; "
##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.17279472 0.015936301 10.842837 3.203038e-23
## homophilyC  -0.03306980 0.024487769 -1.350462 1.779208e-01
## degreeD     -0.04574539 0.006249302 -7.320081 2.450020e-12
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.25 ; "
##                  Estimate  Std. Error     t value     Pr(>|t|)
## (Intercept)  0.1511052406 0.013029224 11.59740896 8.639061e-26
## homophilyC   0.0008594979 0.028204089  0.03047423 9.757099e-01
## degreeD     -0.0477740331 0.007491111 -6.37742981 7.119291e-10
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.5 ; "
##                Estimate Std. Error    t value     Pr(>|t|)
## (Intercept)  0.14045510 0.01490705  9.4220617 1.557726e-18
## homophilyC  -0.02037047 0.02491469 -0.8176088 4.142612e-01
## degreeD     -0.04017664 0.00807308 -4.9766179 1.119228e-06
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.75 ; "
##                Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.14526506 0.016290464  8.9171835 5.566968e-17
## homophilyC   0.01280257 0.023870391  0.5363368 5.921384e-01
## degreeD     -0.04094901 0.006956108 -5.8867701 1.090532e-08
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 1 ; "
##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.14693101 0.014789444  9.934857 3.424056e-20
## homophilyC  -0.03019421 0.021031047 -1.435697 1.521738e-01
## degreeD     -0.03106349 0.005428001 -5.722823 2.626771e-08

## Warning: Removed 17 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 17 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 17 rows containing missing values (`geom_point()`).

## Warning: Removed 17 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 17 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 17 rows containing missing values (`geom_point()`).

## Warning: Removed 20 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 20 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 20 rows containing missing values (`geom_point()`).

## Warning: Removed 23 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 23 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 23 rows containing missing values (`geom_point()`).

## Warning: Removed 14 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 14 rows containing non-finite values (`stat_cor()`).

## Warning: Removed 14 rows containing missing values (`geom_point()`).

## [1] "Regression on proportion of ever-isolated individuals, Control ; "
##                 Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.168216628 0.015666282 10.7374949 9.444767e-23
## homophilyD   0.004659247 0.026161167  0.1780978 8.587749e-01
## degreeD     -0.049197644 0.009203286 -5.3456608 1.871073e-07
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.25 ; "
##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.16838126 0.013416418 12.550388 5.762488e-29
## homophilyD  -0.03165775 0.022777020 -1.389899 1.656637e-01
## degreeD     -0.05109114 0.007758783 -6.584942 2.241807e-10
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.5 ; "
##                Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.14731565 0.015760386  9.3472103 3.128488e-18
## homophilyD  -0.01657153 0.020431273 -0.8110865 4.180124e-01
## degreeD     -0.04372876 0.007501395 -5.8294164 1.544433e-08
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 0.75 ; "
##                Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)  0.15418663 0.016874848  9.1370675 1.475059e-17
## homophilyD   0.01936028 0.019703006  0.9826056 3.266684e-01
## degreeD     -0.04470021 0.007329195 -6.0989247 3.622629e-09
## [1] "Regression on proportion of ever-isolated individuals, degree %ile = 1 ; "
##                 Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept)  0.154872224 0.015556395  9.955535 3.251913e-20
## homophilyD  -0.004801959 0.019178070 -0.250388 8.024690e-01
## degreeD     -0.036807034 0.005913452 -6.224289 1.739021e-09

plot.trends <- 
  data.frame(
    trends.df %>% 
      group_by(round, V, GINI, fractionCoop) %>% 
      summarize_all(list(mean=~mean(., na.rm=TRUE),sd=~sd(., na.rm=TRUE)))
  )

plot.trends$V = factor(plot.trends$V)
plot.trends$GINI = factor(plot.trends$GINI)

for(i in unique(plot.trends$fractionCoop)){
  g.gini = ggplot(data=plot.trends[plot.trends$fractionCoop==i,], aes(x=round,y=gini_mean,group=interaction(GINI,V))) +
  geom_line(aes(color=GINI,linetype=V)) +
  geom_ribbon(aes(ymin = gini_mean - gini_sd, ymax = gini_mean + gini_sd, fill=GINI),alpha=0.3) +
  xlab("Round")+
  ylab("gini") +
  theme_bw() 

  g.gmd = ggplot(data=plot.trends[plot.trends$fractionCoop==i,], aes(x=round,y=gmd_mean,group=interaction(GINI,V))) +
    geom_line(aes(color=GINI,linetype=V)) +
    geom_ribbon(aes(ymin = gmd_mean - gmd_sd, ymax = gmd_mean + gmd_sd, fill=GINI),alpha=0.3) +
    xlab("Round")+
    ylab("gmd") +
    theme_bw() 
  
  g.avg_wealth = ggplot(data=plot.trends[plot.trends$fractionCoop==i,], aes(x=round,y=avg_wealth_mean,group=interaction(GINI,V))) +
    geom_line(aes(color=GINI,linetype=V)) +
    geom_ribbon(aes(ymin = avg_wealth_mean - avg_wealth_sd, ymax = avg_wealth_mean + avg_wealth_sd, fill=GINI),alpha=0.3) +
    xlab("Round")+
    ylab("avg_wealth") +
    theme_bw() 
  
  g.avg_coop = ggplot(data=plot.trends[plot.trends$fractionCoop==i,], aes(x=round,y=avg_coop_mean,group=interaction(GINI,V))) +
    geom_line(aes(color=GINI,linetype=V)) +
    geom_ribbon(aes(ymin = avg_coop_mean - avg_coop_sd, ymax = avg_coop_mean + avg_coop_sd, fill=GINI),alpha=0.3) +
    xlab("Round")+
    ylab("avg_coop") +
    theme_bw() 
  
  g.avg_degree = ggplot(data=plot.trends[plot.trends$fractionCoop==i,], aes(x=round,y=avg_degree_mean,group=interaction(GINI,V))) +
    geom_line(aes(color=GINI,linetype=V)) +
    geom_ribbon(aes(ymin = avg_degree_mean - avg_degree_sd, ymax = avg_degree_mean + avg_degree_sd, fill=GINI),alpha=0.3) +
    xlab("Round")+
    ylab("avg_degree") +
    theme_bw() 
   
  
  
  plot <- ggarrange(g.gini,g.gmd,g.avg_wealth,g.avg_coop,g.avg_degree,common.legend = TRUE,legend="bottom")
  print(annotate_figure(plot, top = text_grob(paste("Degree percentile of nodes assigned to defectors =",i), color = "black", face = "bold", size = 10)))
}

## Warning: Removed 6 rows containing missing values (`geom_line()`).

## Warning: Removed 6 rows containing missing values (`geom_line()`).

## Warning: Removed 6 rows containing missing values (`geom_line()`).

## Warning: Removed 6 rows containing missing values (`geom_line()`).

## Warning: Removed 6 rows containing missing values (`geom_line()`).

fig1 = ggplot(data = df.netIntLowDegree, 
                         aes(x = degreeD, y = homophilyC, color = percentIsolation*100)) + 
  geom_point() +
  scale_x_continuous("Mean degree of defectors") + 
  scale_y_continuous("C-assortativity") +  
  scale_color_viridis(option = "magma") +
  labs(color="Isolated \nindividuals (%)")

fig2 = ggplot(data = df.netIntLowDegree, 
                         aes(x = degreeD, y = homophilyD, color = percentIsolation*100)) + 
  geom_point() +
  scale_x_continuous("Mean degree of defectors") + 
  scale_y_continuous("D-assortativity") +  
  scale_color_viridis(option = "magma") +
  labs(color="Isolated \nindividuals (%)")

fig3 = ggplot(data = df.netIntLowDegree, 
                         aes(x = degreeD, y = heterophily, color = percentIsolation*100)) + 
  geom_point() +
  scale_x_continuous("Mean degree of defectors") + 
  scale_y_continuous("Heterophily") +  
  scale_color_viridis(option = "magma") +
  labs(color="Isolated \nindividuals (%)")

fig4 = ggplot(data = df.netIntLowDegree, 
                         aes(x = degreeC, y = homophilyC, color = percentIsolation*100)) + 
  geom_point() +
  scale_x_continuous("Mean degree of cooperators") + 
  scale_y_continuous("C-assortativity") +  
  scale_color_viridis(option = "magma") +
  labs(color="Isolated \nindividuals (%)")

fig5 = ggplot(data = df.netIntLowDegree, 
                         aes(x = degreeC, y = homophilyD, color = percentIsolation*100)) + 
  geom_point() +
  scale_x_continuous("Mean degree of cooperators") + 
  scale_y_continuous("D-assortativity") +  
  scale_color_viridis(option = "magma") +
  labs(color="Isolated \nindividuals (%)")


fig6 = ggplot(data = df.netIntLowDegree, 
                         aes(x = degreeC, y = heterophily, color = percentIsolation*100)) + 
  geom_point() +
  scale_x_continuous("Mean degree of cooperators") + 
  scale_y_continuous("Heterophily") +  
  scale_color_viridis(option = "magma") +
  labs(color="Isolated \nindividuals (%)")

fig7 = ggplot(data = df.netIntLowDegree, 
                         aes(x = degreeC, y = degreeD, color = percentIsolation*100)) + 
  geom_point() +
  scale_x_continuous("Mean degree of cooperators") + 
  scale_y_continuous("Mean degree of defectors") +  
  scale_color_viridis(option = "magma") +
  labs(color="Isolated \nindividuals (%)")

print(ggarrange(fig1,fig2,fig3,fig4,fig5,fig6,fig7,common.legend = TRUE,legend="right"))

## Warning: Removed 43 rows containing missing values (`geom_point()`).
## Removed 43 rows containing missing values (`geom_point()`).

## Warning: Removed 91 rows containing missing values (`geom_point()`).

## Warning: Removed 12 rows containing missing values (`geom_point()`).

## Warning: Removed 33 rows containing missing values (`geom_point()`).

## Warning: Removed 93 rows containing missing values (`geom_point()`).

## Warning: Removed 12 rows containing missing values (`geom_point()`).
## Removed 12 rows containing missing values (`geom_point()`).

fig1 = ggplot(data = df.netIntLowDegree, 
                         aes(x = degreeC, y = percentIsolation*100)) + 
  geom_point() +
  scale_x_continuous("Mean degree of cooperators") + 
  scale_y_continuous("Isolated individuals (%)") 


fig2 = ggplot(data = df.netIntLowDegree, 
                         aes(x = degreeD, y = percentIsolation*100)) + 
  geom_point() +
  scale_x_continuous("Mean degree of defectors") + 
  scale_y_continuous("Isolated individuals (%)") 


fig3 = ggplot(data = df.netIntLowDegree, 
                         aes(x = homophilyC, y = percentIsolation*100)) + 
  geom_point() +
  scale_x_continuous("C-assortativity") + 
  scale_y_continuous("Isolated individuals (%)") 



fig4 = ggplot(data = df.netIntLowDegree, 
                         aes(x = homophilyD, y = percentIsolation*100)) + 
  geom_point() +
  scale_x_continuous("D-assortativity") + 
  scale_y_continuous("Isolated individuals (%)") 



print(ggarrange(fig1,fig2,fig3,fig4,common.legend = TRUE,legend="right"))

## Warning: Removed 2 rows containing missing values (`geom_point()`).

## Warning: Removed 10 rows containing missing values (`geom_point()`).

## Warning: Removed 33 rows containing missing values (`geom_point()`).

## Warning: Removed 91 rows containing missing values (`geom_point()`).

## Warning: Removed 2 rows containing missing values (`geom_point()`).

## Warning: Removed 10 rows containing missing values (`geom_point()`).

## Warning: Removed 33 rows containing missing values (`geom_point()`).

## Warning: Removed 91 rows containing missing values (`geom_point()`).

reg.isolation = glm(percentIsolation*100 ~ degreeC + degreeD + homophilyC + homophilyD + heterophily, data=df.netIntLowDegree, family = gaussian(link = "identity"))
summary(reg.isolation)

## 
## Call:
## glm(formula = percentIsolation * 100 ~ degreeC + degreeD + homophilyC + 
##     homophilyD + heterophily, family = gaussian(link = "identity"), 
##     data = df.netIntLowDegree)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -15.740   -6.261   -2.530    4.719   38.258  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  18.3649     0.8200  22.397  < 2e-16 ***
## degreeC      -2.1114     0.4882  -4.325 1.64e-05 ***
## degreeD      -3.5120     0.4013  -8.751  < 2e-16 ***
## homophilyC    0.7154     1.3456   0.532    0.595    
## homophilyD   -0.1944     1.0394  -0.187    0.852    
## heterophily  -0.9728     3.3468  -0.291    0.771    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 73.58608)
## 
##     Null deviance: 119633  on 1375  degrees of freedom
## Residual deviance: 100813  on 1370  degrees of freedom
##   (124 observations deleted due to missingness)
## AIC: 9827.6
## 
## Number of Fisher Scoring iterations: 2

#variance inflation factor
car::vif(reg.isolation)

##     degreeC     degreeD  homophilyC  homophilyD heterophily 
##    2.475714    2.104270    1.649504    1.398580    2.958421

reg.isolation = glm(percentIsolation*100 ~ degreeD + homophilyD, data=df.netIntLowDegree, family = gaussian(link = "identity"))
summary(reg.isolation)

## 
## Call:
## glm(formula = percentIsolation * 100 ~ degreeD + homophilyD, 
##     family = gaussian(link = "identity"), data = df.netIntLowDegree)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -14.991   -6.527   -2.779    4.837   39.840  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  14.9908     0.6318  23.727   <2e-16 ***
## degreeD      -3.9419     0.2920 -13.501   <2e-16 ***
## homophilyD   -1.1593     0.9332  -1.242    0.214    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 75.63175)
## 
##     Null deviance: 122497  on 1408  degrees of freedom
## Residual deviance: 106338  on 1406  degrees of freedom
##   (91 observations deleted due to missingness)
## AIC: 10099
## 
## Number of Fisher Scoring iterations: 2

#variance inflation factor
car::vif(reg.isolation)

##    degreeD homophilyD 
##   1.101326   1.101326

#double machine learning
library(DoubleML)
library(mlr3)
library(mlr3learners)
set.seed(3141)


##degreeC
dml_data = DoubleMLData$new(df.netIntLowDegree[complete.cases(df.netIntLowDegree[c("percentIsolation","degreeC","degreeD","homophilyC","homophilyD","heterophily")]),],
                             y_col = "percentIsolation",
                             d_cols = "degreeC",
                             x_cols = c("degreeD","homophilyC","homophilyD","heterophily"))
print(dml_data)

## ================= DoubleMLData Object ==================
## 
## 
## ------------------ Data summary      ------------------
## Outcome variable: percentIsolation
## Treatment variable(s): degreeC
## Covariates: degreeD, homophilyC, homophilyD, heterophily
## Instrument(s): 
## No. Observations: 1376

# surpress messages from mlr3 package during fitting
lgr::get_logger("mlr3")$set_threshold("warn")

learner = lrn("regr.ranger", num.trees=500, mtry=floor(sqrt(4)), max.depth=5, min.node.size=2)
ml_l = learner$clone()
ml_m = learner$clone()

obj_dml_plr = DoubleMLPLR$new(dml_data, ml_l=ml_l, ml_m=ml_m)
obj_dml_plr$fit()
print(obj_dml_plr)

## ================= DoubleMLPLR Object ==================
## 
## 
## ------------------ Data summary      ------------------
## Outcome variable: percentIsolation
## Treatment variable(s): degreeC
## Covariates: degreeD, homophilyC, homophilyD, heterophily
## Instrument(s): 
## No. Observations: 1376
## 
## ------------------ Score & algorithm ------------------
## Score function: partialling out
## DML algorithm: dml2
## 
## ------------------ Machine learner   ------------------
## ml_l: regr.ranger
## ml_m: regr.ranger
## 
## ------------------ Resampling        ------------------
## No. folds: 5
## No. repeated sample splits: 1
## Apply cross-fitting: TRUE
## 
## ------------------ Fit summary       ------------------
##  Estimates and significance testing of the effect of target variables
##         Estimate. Std. Error t value Pr(>|t|)    
## degreeC -0.020275   0.004813  -4.213 2.52e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##degreeD
dml_data = DoubleMLData$new(df.netIntLowDegree[complete.cases(df.netIntLowDegree[c("percentIsolation","degreeC","degreeD","homophilyC","homophilyD","heterophily")]),],
                             y_col = "percentIsolation",
                             d_cols = "degreeD",
                             x_cols = c("degreeC","homophilyC","homophilyD","heterophily"))
print(dml_data)

## ================= DoubleMLData Object ==================
## 
## 
## ------------------ Data summary      ------------------
## Outcome variable: percentIsolation
## Treatment variable(s): degreeD
## Covariates: degreeC, homophilyC, homophilyD, heterophily
## Instrument(s): 
## No. Observations: 1376

# surpress messages from mlr3 package during fitting
lgr::get_logger("mlr3")$set_threshold("warn")

learner = lrn("regr.ranger", num.trees=500, mtry=floor(sqrt(4)), max.depth=5, min.node.size=2)
ml_l = learner$clone()
ml_m = learner$clone()

obj_dml_plr = DoubleMLPLR$new(dml_data, ml_l=ml_l, ml_m=ml_m)
obj_dml_plr$fit()
print(obj_dml_plr)

## ================= DoubleMLPLR Object ==================
## 
## 
## ------------------ Data summary      ------------------
## Outcome variable: percentIsolation
## Treatment variable(s): degreeD
## Covariates: degreeC, homophilyC, homophilyD, heterophily
## Instrument(s): 
## No. Observations: 1376
## 
## ------------------ Score & algorithm ------------------
## Score function: partialling out
## DML algorithm: dml2
## 
## ------------------ Machine learner   ------------------
## ml_l: regr.ranger
## ml_m: regr.ranger
## 
## ------------------ Resampling        ------------------
## No. folds: 5
## No. repeated sample splits: 1
## Apply cross-fitting: TRUE
## 
## ------------------ Fit summary       ------------------
##  Estimates and significance testing of the effect of target variables
##         Estimate. Std. Error t value Pr(>|t|)    
## degreeD -0.030903   0.004077   -7.58 3.46e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##homophilyC
dml_data = DoubleMLData$new(df.netIntLowDegree[complete.cases(df.netIntLowDegree[c("percentIsolation","degreeC","degreeD","homophilyC","homophilyD","heterophily")]),],
                             y_col = "percentIsolation",
                             d_cols = "homophilyC",
                             x_cols = c("degreeC","degreeD","homophilyD","heterophily"))
print(dml_data)

## ================= DoubleMLData Object ==================
## 
## 
## ------------------ Data summary      ------------------
## Outcome variable: percentIsolation
## Treatment variable(s): homophilyC
## Covariates: degreeC, degreeD, homophilyD, heterophily
## Instrument(s): 
## No. Observations: 1376

# surpress messages from mlr3 package during fitting
lgr::get_logger("mlr3")$set_threshold("warn")

learner = lrn("regr.ranger", num.trees=500, mtry=floor(sqrt(4)), max.depth=5, min.node.size=2)
ml_l = learner$clone()
ml_m = learner$clone()

obj_dml_plr = DoubleMLPLR$new(dml_data, ml_l=ml_l, ml_m=ml_m)
obj_dml_plr$fit()
print(obj_dml_plr)

## ================= DoubleMLPLR Object ==================
## 
## 
## ------------------ Data summary      ------------------
## Outcome variable: percentIsolation
## Treatment variable(s): homophilyC
## Covariates: degreeC, degreeD, homophilyD, heterophily
## Instrument(s): 
## No. Observations: 1376
## 
## ------------------ Score & algorithm ------------------
## Score function: partialling out
## DML algorithm: dml2
## 
## ------------------ Machine learner   ------------------
## ml_l: regr.ranger
## ml_m: regr.ranger
## 
## ------------------ Resampling        ------------------
## No. folds: 5
## No. repeated sample splits: 1
## Apply cross-fitting: TRUE
## 
## ------------------ Fit summary       ------------------
##  Estimates and significance testing of the effect of target variables
##            Estimate. Std. Error t value Pr(>|t|)
## homophilyC -0.002112   0.012403   -0.17    0.865

##homophilyD
dml_data = DoubleMLData$new(df.netIntLowDegree[complete.cases(df.netIntLowDegree[c("percentIsolation","degreeC","degreeD","homophilyC","homophilyD","heterophily")]),],
                             y_col = "percentIsolation",
                             d_cols = "homophilyD",
                             x_cols = c("degreeC","degreeD","homophilyC","heterophily"))
print(dml_data)

## ================= DoubleMLData Object ==================
## 
## 
## ------------------ Data summary      ------------------
## Outcome variable: percentIsolation
## Treatment variable(s): homophilyD
## Covariates: degreeC, degreeD, homophilyC, heterophily
## Instrument(s): 
## No. Observations: 1376

# surpress messages from mlr3 package during fitting
lgr::get_logger("mlr3")$set_threshold("warn")

learner = lrn("regr.ranger", num.trees=500, mtry=floor(sqrt(4)), max.depth=5, min.node.size=2)
ml_l = learner$clone()
ml_m = learner$clone()

obj_dml_plr = DoubleMLPLR$new(dml_data, ml_l=ml_l, ml_m=ml_m)
obj_dml_plr$fit()
print(obj_dml_plr)

## ================= DoubleMLPLR Object ==================
## 
## 
## ------------------ Data summary      ------------------
## Outcome variable: percentIsolation
## Treatment variable(s): homophilyD
## Covariates: degreeC, degreeD, homophilyC, heterophily
## Instrument(s): 
## No. Observations: 1376
## 
## ------------------ Score & algorithm ------------------
## Score function: partialling out
## DML algorithm: dml2
## 
## ------------------ Machine learner   ------------------
## ml_l: regr.ranger
## ml_m: regr.ranger
## 
## ------------------ Resampling        ------------------
## No. folds: 5
## No. repeated sample splits: 1
## Apply cross-fitting: TRUE
## 
## ------------------ Fit summary       ------------------
##  Estimates and significance testing of the effect of target variables
##            Estimate. Std. Error t value Pr(>|t|)
## homophilyD -0.004637   0.008776  -0.528    0.597

##heterophily
dml_data = DoubleMLData$new(df.netIntLowDegree[complete.cases(df.netIntLowDegree[c("percentIsolation","degreeC","degreeD","homophilyC","homophilyD","heterophily")]),],
                             y_col = "percentIsolation",
                             d_cols = "heterophily",
                             x_cols = c("degreeC","degreeD","homophilyC","homophilyD"))
print(dml_data)

## ================= DoubleMLData Object ==================
## 
## 
## ------------------ Data summary      ------------------
## Outcome variable: percentIsolation
## Treatment variable(s): heterophily
## Covariates: degreeC, degreeD, homophilyC, homophilyD
## Instrument(s): 
## No. Observations: 1376

# surpress messages from mlr3 package during fitting
lgr::get_logger("mlr3")$set_threshold("warn")

learner = lrn("regr.ranger", num.trees=500, mtry=floor(sqrt(4)), max.depth=5, min.node.size=2)
ml_l = learner$clone()
ml_m = learner$clone()

obj_dml_plr = DoubleMLPLR$new(dml_data, ml_l=ml_l, ml_m=ml_m)
obj_dml_plr$fit()
print(obj_dml_plr)

## ================= DoubleMLPLR Object ==================
## 
## 
## ------------------ Data summary      ------------------
## Outcome variable: percentIsolation
## Treatment variable(s): heterophily
## Covariates: degreeC, degreeD, homophilyC, homophilyD
## Instrument(s): 
## No. Observations: 1376
## 
## ------------------ Score & algorithm ------------------
## Score function: partialling out
## DML algorithm: dml2
## 
## ------------------ Machine learner   ------------------
## ml_l: regr.ranger
## ml_m: regr.ranger
## 
## ------------------ Resampling        ------------------
## No. folds: 5
## No. repeated sample splits: 1
## Apply cross-fitting: TRUE
## 
## ------------------ Fit summary       ------------------
##  Estimates and significance testing of the effect of target variables
##             Estimate. Std. Error t value Pr(>|t|)
## heterophily  -0.02018    0.03342  -0.604    0.546

fig1 = ggplot(data = df.netIntLowDegree, 
                         aes(x = degreeD, y = homophilyC, color = avgCoopFinal*100)) + 
  geom_point() +
  scale_x_continuous("Mean degree of defectors") + 
  scale_y_continuous("C-assortativity") +  
  scale_color_viridis(option = "magma") +
  labs(color="Cooperation in \nfinal round (%)")

fig2 = ggplot(data = df.netIntLowDegree, 
                         aes(x = degreeD, y = homophilyD, color = avgCoopFinal*100)) + 
  geom_point() +
  scale_x_continuous("Mean degree of defectors") + 
  scale_y_continuous("D-assortativity") +  
  scale_color_viridis(option = "magma") +
  labs(color="Cooperation in \nfinal round (%)")

fig3 = ggplot(data = df.netIntLowDegree, 
                         aes(x = degreeD, y = heterophily, color = avgCoopFinal*100)) + 
  geom_point() +
  scale_x_continuous("Mean degree of defectors") + 
  scale_y_continuous("Heterophily") +  
  scale_color_viridis(option = "magma") +
  labs(color="Cooperation in \nfinal round (%)")

fig4 = ggplot(data = df.netIntLowDegree, 
                         aes(x = degreeC, y = homophilyC, color = avgCoopFinal*100)) + 
  geom_point() +
  scale_x_continuous("Mean degree of cooperators") + 
  scale_y_continuous("C-assortativity") +  
  scale_color_viridis(option = "magma") +
  labs(color="Cooperation in \nfinal round (%)")

fig5 = ggplot(data = df.netIntLowDegree, 
                         aes(x = degreeC, y = homophilyD, color = avgCoopFinal*100)) + 
  geom_point() +
  scale_x_continuous("Mean degree of cooperators") + 
  scale_y_continuous("D-assortativity") +  
  scale_color_viridis(option = "magma") +
  labs(color="Cooperation in \nfinal round (%)")


fig6 = ggplot(data = df.netIntLowDegree, 
                         aes(x = degreeC, y = heterophily, color = avgCoopFinal*100)) + 
  geom_point() +
  scale_x_continuous("Mean degree of cooperators") + 
  scale_y_continuous("Heterophily") +  
  scale_color_viridis(option = "magma") +
  labs(color="Cooperation in \nfinal round (%)")

fig7 = ggplot(data = df.netIntLowDegree, 
                         aes(x = degreeC, y = degreeD, color = avgCoopFinal*100)) + 
  geom_point() +
  scale_x_continuous("Mean degree of cooperators") + 
  scale_y_continuous("Mean degree of defectors") +  
  scale_color_viridis(option = "magma") +
  labs(color="Cooperation in \nfinal round (%)")

print(ggarrange(fig1,fig2,fig3,fig4,fig5,fig6,fig7,common.legend = TRUE,legend="right"))

## Warning: Removed 43 rows containing missing values (`geom_point()`).
## Removed 43 rows containing missing values (`geom_point()`).

## Warning: Removed 91 rows containing missing values (`geom_point()`).

## Warning: Removed 12 rows containing missing values (`geom_point()`).

## Warning: Removed 33 rows containing missing values (`geom_point()`).

## Warning: Removed 93 rows containing missing values (`geom_point()`).

## Warning: Removed 12 rows containing missing values (`geom_point()`).
## Removed 12 rows containing missing values (`geom_point()`).

Several explanations

More beneficial to break links

Mathematically, the expected benefit of breaking/not making new links only depends on the total number of cooperators and defectors, not the number of C and D around an ego

Altruism: making/retaining links even when there is no benefit

Degree changes in initial defectors are no different across conditions

Evolution of cooperation

This seems to occur regardless of initial defector/cooperator position (the same number of C's become D's, and the same number of D's become C's)

D has more links to start

This seems to be the most plausible explanation. Regardless of the network position, initial D's ends up becoming C's in most cases, but not without a cost: many D's tend to lose links first before reaching equilibrium, so networks that have the greatest D degrees may have the least isolation.